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Preface 


This  report  describes  the  evaluation  and  application  of  the  Department 
of  Defense  Groundwater  Modeling  System  (GMS,  Version  1.2)  at 
Columbus  Air  Force  Base,  Columbus,  MS.  The  project  was  conducted  as 
part  of  the  U.  S.  Army  Engineer  Waterways  Experiment  Station  (WES) 
groundwater  model  development  and  application  program.  The  primary 
objective  of  this  work  was  to  evaluate  GMS/FEMWATER  and 
GMS/MODFLOW-MT3D  using  data  from  the  MADE-2  groundwater  flow 
and  transport  experiment.  The  project  was  performed  as  a  component  of 
the  WES  Groundwater  Modeling  Program  and  was  funded  by  the  U.S. 
Army  Environmental  Center  (AEC).  Mr.  Ira  May  is  the  AEC  technical 
monitor  for  the  project.  Field  data  were  provided  by  Dr.  Thomas  Stauffer, 
Armstrong  Laboratory  Environics,  Tyndall  Air  Force  Base,  FL. 

The  study  was  conducted  under  the  direct  supervision  of  Dr.  Mark  S. 
Dortch,  Chief,  Water  Quality  and  Contaminant  Modeling  Branch 
(WQCMB),  and  under  the  general  supervision  of  Dr.  Richard  E.  Price, 
Chief,  Environmental  Processes  and  Effects  Division,  and  Dr.  John 
Harrison,  Director,  Environmental  Laboratory.  This  report  was  written  by 
Dr.  Mansour  Zakikhani,  WQCMB,  Mr.  Kenneth  C.  Black,  Environmental 
Consulting  Engineers,  and  Mr.  Jeffrey  A.  Gerald,  AScI  Corporation. 

Acknowledgment  is  made  to  Mr.  Chris  McGrath,  WQCMB,  and 
Dr.  Fred  Tracy,  WES  Information  Technology  Laboratory,  for  their  review 
and  valuable  suggestions. 

The  work  was  coordinated  by  Dr.  Jeffery  P.  Holland,  Director, 
Computational  Hydraulics  Institute,  WES  Coastal  and  Hydraulics 
Laboratory,  and  WES  Groundwater  Modeling  Program  Manager. 

At  the  time  of  publication  of  this  report,  Dr.  Robert  W.  Whalin  was 
Director  of  WES. 
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1  Introduction 


The  Department  of  Defense  Groundwater  Modeling  System  (GMS)  is  a 
comprehensive  computer  graphical  system.  The  GMS  has  been  developed 
as  part  of  an  ongoing  multi-disciplinary  research  effort.  The  purpose  of  the 
research  initiative  was  to  develop  a  user  environment  where  existing 
groundwater  models  can  be  used  efficiently.  The  GMS  includes  modeling 
tools  to  facilitate  site  characterization,  model  conceptualization,  mesh  and 
grid  generation,  geostatistical  computations,  and  post-processing.  Version 
1.2  of  the  GMS  provides  an  interface  for  subsurface  models  such  as 
FEMWATER  and  MODFLOW-MT3D.  As  part  of  a  project  on  the  GMS 
evaluation,  development,  and  improvement,  the  GMS  has  been  applied  to 
selected  field  problems.  This  rigorous  evaluation  of  GMS  capabilities  is 
intended  to  identify  additional  numerical  tools  needed  for  advancement  of 
the  GMS. 

One  of  the  sites  which  has  received  significant  attention  and  provides 
extensive  data  sets  is  located  on  the  Columbus  Air  Force  Base  (CAFB), 
Mississippi.  CAFB  site  is  the  location  of  two  macrodispersion  experiments 
(MADE-1  and  MADE-2)  and  the  present  ongoing  natural  attenuation  study 
(NATS).  MADE-1  involved  the  injection  and  monitoring  of  bromide  and 
three  flourobenzoate  tracers  between  October  1986  and  1988  (Boggs  et  al. 
1988;  Boggs  et  al.  1990).  MADE-2  (Boggs  et  al.  1993)  studied  physical, 
chemical,  and  biological  processes  affecting  transport  of  two  conservative 
chemicals  and  four  reactive  hydrocarbons  injected  as  solutes.  The 
flowmeter-derived  hydraulic  conductivity  profile  of  the  site  illustrated  the 
extreme  heterogeneity  of  the  aquifer.  The  hydraulic  head  field  at  the  site  has 
exhibited  complex  temporal  and  spatial  variability  produced  by  the  large 
heterogeneity  of  the  aquifer  and  large  seasonal  fluctuations  of  the  water 
table.  The  MADE-2  study  was  initiated  with  a  2-day  injection  of  the 
dissolved  chemicals  beginning  June  26,  1990.  Chemical  concentration 
distributions  were  monitored  using  an  extensive  network  of  saturated  zone 
multilevel  samples  at  1-  to  3-month  intervals  over  a  period  of  15  months. 
Detailed  information  on  field  monitoring  may  be  found  in  Boggs  et  al. 
(1993).  Although  extensive  data  collections  and  field  monitoring  have  been 
conducted  at  the  site,  few  numerical  modeling  studies  have  been  reported 
(Gray  and  Dale  1995).  As  part  of  a  groundwater  model  development  at  the 
U.  S.  Army  Engineer  Waterways  Experiment  Station  selected  numerical 
models  in  the  GMS  including  FEMWATER  and  MODFLOW-MT3D  have 
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been  applied  to  the  site  using  MADE-2  data.  In  this  report,  the  capability  of 
the  GMS  Version  1.2  (Brigham  Young  University  1995)  will  be  discussed 
and  illustrated  by  applying  it  to  this  complex  heterogeneous  aquifer.  The 
report  covers  two  different  numerical  modeling  approaches:  finite  element 
(FEMWATER)  and  finite  difference  (MODFLOW-MT3D). 


Purpose 

The  primary  purpose  of  this  demonstration  project  was  to  determine  the 
current  strengths  and  weaknesses  of  the  GMS  (Version  1.2)  and  its 
groundwater  models  through  a  field  modeling  application.  Although  an 
attempt  was  made  to  calibrate  the  flow  and  transport  models  to  measured 
data,  the  objective  of  this  project  was  not  to  create  a  highly  refined  model 
of  the  MADE-2  experiment.  The  MADE-2  tracer  experiment  was  selected 
for  the  application  because  the  test  is  well-documented  and  the 
hydrogeologic  database  of  the  experiment  was  available.  Additionally, 
accurate  numerical  simulation  of  this  experiment  is  generally  considered  to 
be  a  difficult  task.1  For  this  reason,  this  experiment  should  provide  ample 
opportunity  to  test  many  of  the  features  of  the  GMS/FEMWATER  and 
GMS/MODFLOW-  MT3D  software  system. 


Scope  of  Report 

This  report  is  divided  into  two  main  parts.  In  the  first  part  of  the  report 
(Chapters  2,  3,  and  4)  the  reader  can  find  information  on  general 
characteristics  of  the  site,  and  the  application,  modification,  and  merits  of 
the  GMS/FEMWATER  software  system.  The  hydrogeologic  setting  of  the 
MADE-2  experiment  and  the  hydrogeologic  database  are  described.  Results 
of  flow  and  transport  simulations  are  presented  followed  by  conclusions 
determined  from  the  modeling  exercise.  In  the  second  part  of  the  report 
(Chapters  5  and  6),  the  application  of  GMS/MODFLOW-MT3D  to 
MADE-2  sites  is  discussed  and  illustrated.  The  modeling  approaches  used 
to  create  the  flow  and  transport  models  are  discussed  in  detail.  Interpolated 
hydraulic  conductivity  values  for  each  layer  of  the  model  grid  are  presented 
in  Appendix  A.  Appendix  B  contains  average  head  contours  for  both  the 
shallow  and  deep  parts  of  the  aquifer  for  each  day  of  the  available  data. 
Gauge  plots  of  transient  flow  simulations  are  presented  in  Appendix  C. 


1  Personal  Communication  on  the  history  of  work  conducted  at  the  MADE  site,  1996,  Steve 
Young,  Environmental  Consulting  Engineers,  Inc.,  Knoxville,  TN. 
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2  Overview  of  Site 
Hydrogeology 


Site  Description 

The  study  area  is  located  on  the  Columbus  Air  Force  Base  in  Columbus, 
MS  (Figure  1).  The  site  is  known  as  MADE  because  of  a  macro  dispersion 
experiment  conducted  at  the  site.  The  air  base  is  situated  above  a  shallow 
unconfined  aquifer,  which  consists  of  fluvial  deposits  underlain  by  a  Eutaw 
Formation  marine  sediment  aquitard  composed  of  clays,  silts,  and 
fine-grained  sands.  The  average  thickness  of  the  aquifer  is  approximately 
11m  and  it  is  composed  of  poorly  sorted  to  well-sorted  sandy  gravel  and 
gravelly  sand  with  some  silt  and  clay.  The  aquifer  contains  several  clay 
lenses  and  there  is  apparent  evidence  of  a  prior  meandering  channel 
crossing  the  deposit  (Figure  2).  The  surface  of  the  site  is  gently  sloping 
northward  from  an  elevation  in  the  southeastern  corner  of  66.5  m  above 
mean  sea  level  (A.MSL)  to  64.6  m  (A.MSL)  in  the  northern  corner.  Soil 
facies,  determined  from  a  limited  terrace  exposure,  indicate  irregular  lenses 
and  layers  with  horizontal  dimensions  of  8  m  and  vertical  dimensions  of 
less  than  1  m.  The  piezometric  head  and  hydraulic  gradients  at  the  site 
exhibit  complex  spatial  and  temporal  variability.  The  saturated  thickness  of 
the  aquifer  varies  20  to  30  percent  between  the  low  and  high  water  table 
seasons.  The  area  near  the  injection  wells  has  steep  horizontal  gradients 
while  the  area  furthest  north  contains  less  of  a  horizontal  gradient  but  more 
of  a  vertical  gradient. 

The  mean  horizontal  hydraulic  conductivity  near  the  injection  site 
(southernmost  end)  and  at  the  end  (northernmost)  of  the  test  area  is 
approximately  0.86  m/day  with  significant  variation  between  these  two 
areas  (e.g.  local  areas  in  excess  of  735  m/day  are  indicated).  The  central 
region  with  relatively  high  mean  hydraulic  conductivities  generally 
corresponds  to  the  former  river  meander  channel  that  trends  southwest  to 
northeast  and  lies  at  an  approximate  45-deg  angle  from  the  direction  of 
groundwater  flow.  Vertical  variability  of  hydraulic  conductivity  within  a 
measurement  well  typically  varies  by  two  to  four  orders  of  magnitude. 
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Figure  1. 
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SCALE  KILOMETERS 

Site  location  map  with  regional  features  such  as  rivers  and  oxbow 
channels 
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Concave  Paleochanne! 


Figure  2.  Aerial  photo  showing  the  MADE-2  injection  site  and  the  location  of  the 
meander  channel,  with  an  orthographic  projection  showing  the  fluvial 
deposits  of  a  typical  channel 
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Analysis  of  the  soil  mineralogy  for  the  fine-grained  deposits  of  the  site 
indicate  quartz,  muscovite,  potassium  feldspar,  and  clay  minerals 
(predominantly  kaolinites  and  illites,  with  minor  amounts  of  vermiculite 
and  montmorillonite).  The  sand  and  gravel  are  predominantly  composed  of 
feldspar,  mica,  and  quartz. 

The  climate  at  the  MADE  site  is  temperate  with  a  mean  annual  air 
temperature  of  63°F.  Most  precipitation  at  the  site  occurs  during  the  early 
spring  and  winter,  with  the  least  amount  during  the  early  fall  and  summer. 

A  more  detailed  description  of  the  site  can  be  found  in  several  references 
from  previous  studies  (e.g.  Rehfelt  et  al.  1992,  Boggs  et  al.  1990). 


Description  of  the  Experiment 

There  are  numerous  publications  that  describe  the  site  conditions  that 
were  measured  either  during  the  MADE-1  tracer  experiment  (  Boggs  et  al. 
1988,  Boggs  and  Young  1988)  or  the  MADE-2  tracer  experiment  (Electric 
Power  Research  Institute  (EPRI)  1993a).  The  MADE-1  experiment, 
conducted  between  October  1986  and  June  1988,  examined  the 
scale-dependence  of  the  macrodispersion  coefficients  for  nonreactive 
tracers  (EPRI  1991).  Bromide  was  the  primary  tracer  used  during  this 
experiment.  The  MADE-2  experiment,  conducted  between  June  1990  and 
September  1991,  included  nonreactive,  reactive,  and  degrading  chemicals. 
Four  organic  chemicals,  including  14C-labeled  p-xylene,  o-dichlorobenzene, 
naphthalene,  and  benzene  were  used  along  with  the  conservative  tracer 
tritium.  FEMWATER  was  used  to  simulate  the  transport  of  the  tritium 
tracer.  MODFLOW/MT3D  was  used  to  simulate  the  transport  of  all  the 
chemicals  mentioned  above. 

During  June  27-28,  1990,  9.7  m3  of  the  solution  was  injected  through 
five  wells  spaced  1  m  apart  having  0.6-m  screen  sections.  The  vertical 
injection  interval  ranged  from  an  elevation  of  57.5  m  to  58.1  m  above  mean 
sea  level.  The  initial  mean  aqueous  concentrations  are  55,610  pCi/mL  for 
tritium,  2,770  pCi/ml  for  14C-xylene,  68.1  mg/L  for  benzene,  51.5  mg/L  for 
p-xylene,  7.23  mg/L  for  naphthalene,  and  32.8  mg/L  for  o-dichlorobenzene 
(EPRI  1993b).  Monthly  water  levels  were  recorded  during  the  experiment  at 
numerous  wells,  and  water  levels  at  twelve  wells  were  continuously 
monitored  throughout  the  experiment.  Tracer  concentrations  were  sampled 
throughout  the  experiment.  The  available  complete  concentration  data 
correspond  to  27,  132,  224,  and  328  days  after  injection. 

Some  data  from  the  MADE-1  experiment  were  used  in  this  work  (e.g., 
borehole  flowmeter  results,  well  coordinates).  One  difference  between  the 
experiments  is  that  a  coordinate  transformation  is  required  to  convert  the 
MADE-1  coordinates  to  MADE-2  coordinates  (EPRI  1993b,  page  3-1).  The 
MADE-2  coordinate  system  was  used  for  the  FEMWATER  applications. 
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Hydrogeologic  and  experimental  data 

The  hydrogeologic  and  chemical  data  presented  in  this  report  have  been 
acquired  from  numerous  published  reports  and  the  MADE-2  database  (EPRI 
1993b).  Descriptions  of  site  conditions  were  obtained  from  several  reports 
(EPRI  1989,  1990,  1993a, b).  Details  of  the  MADE-2  tracer  test  were 
obtained  from  several  sources  (EPRI  1993). 1  The  digital  data,  such  as  the 
borehole  flowmeter  hydraulic  conductivities,  have  been  used  within  GMS  to 
create  the  data  needed  for  FEMWATER.  Digital  data  utilized  for  this  model 
include  but  are  not  limited  to: 

a .  Borehole  flowmeter  hydraulic  conductivities. 

b .  Water  level  data  collected  on  a  monthly  basis. 

c.  Continuously  recorded  groundwater  data. 

d .  Tracer  data  collected  at  all  locations  at  the  selected  times. 


Hydraulic  conductivity 

A  total  of  2,018  individual  hydraulic  conductivities  were  available  in  the 
data  set  (Rehfeldt  et  al.  1989a,  1989b).  These  data  were  analyzed  and 
converted  to  the  form  that  could  be  used  for  GMS/FEMWATER  simulations. 
How  these  data  were  incorporated  into  the  model  is  discussed  in  Chapter  3. 


Monthly  water  levels 

Monthly  water  levels  were  imported  to  a  Microsoft  Excel  spreadsheet 
and  used  to  analyze  trends  in  the  data.  Horizontal  and  vertical  hydraulic 
gradients,  transient  water  level  responses,  and  snap-shot  water  level  rounds 
were  studied  using  this  data  set. 


Continuously  recorded  groundwater  levels 

Water  levels  from  16  wells  were  continuously  recorded  for  449  days 
during  the  experiment.  This  information  was  analyzed  to  determine  the 
impact  of  transient  water  level  changes  on  tracer  migration.  The  water 
levels  collected  came  from  eight  vertically  nested  well  pairs  including: 
P-27,  P-44,  P-53,  P-54,  P-55,  P-60,  P-61,  and  P-8  (Figure  3  and  Table  1). 
The  “A”  series  wells  are  the  shallow  wells  and  the  “B”  series  wells  are  the 
deep  wells. 


1  Personal  communication  on  the  history  of  work  conducted  at  the  MADE  site,  1996,  Steve 
Young,  Environmental  Consulting  Engineers,  Inc.,  Knoxville,  TN. 
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Figure  3.  Well  locations  with  nested  wells  labeled  “A/B” 
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Table  1 

Monitoring  Wells 

ID 

X  (m) 

Y  (m) 

Well  ID 

ID 

X(m) 

Y  (m) 

Well  ID 

1 

299.91 

141.05 

P-2 

19 

3.16 

147.63 

P-27A/B 

2 

-86.10 

2.56 

P-3 

20 

32.44 

-25.81 

P-34A/B 

3 

33.56 

-92.08 

P-16 

21 

201.85 

171.53 

P-36A/B/C/D 

4 

32.61 

-157.50 

P-20 

22 

5.55 

71.69 

P-42A/B 

5 

-34.54 

-124.09 

P-23 

23 

43.33 

94.53 

P-43A/B 

6 

130.21 

-345.63 

P-28 

24 

-59.80 

67.92 

P-44A/B 

7 

242.34 

-66.98 

P-29 

25 

52.20 

47.77 

P-46A/B 

8 

-11.41 

83.28 

P-40 

26 

-4.92 

-10.33 

P-53A/B 

9 

1.62 

56.17 

P-41 

27 

-20.13 

35.08 

P-54A/B 

10 

-30.84 

10.38 

P-45 

28 

29.26 

22.95 

P-55A/B 

11 

-85.85 

190.51 

P-52 

29 

-32.26 

113.18 

P-56A/B 

12 

217.69 

366.34 

P-5R 

30 

29.24 

127.12 

P-57A/B 

13 

95.13 

123.75 

P-8A/B 

31 

-69.15 

139.66 

P-58A/B 

14 

112.22 

-201.80 

P-1 

32 

-25.70 

179.32 

P-59A/B 

15 

103.20 

5.73 

P-10A/B 

33 

29.34 

225.10 

P-60A/B 

16 

-19.40 

-50.16 

P-22A/B 

34 

-47.80 

263.96 

P-61  A/B 

17 

22.61 

6.42 

P-24A/B 

35 

51.93 

173.77 

P-62A/B 

18 

-23.30 

66.57 

P-25A/B/C 

Tracer  data 

All  tracer  concentration  data  were  imported  into  a  Microsoft  Excel 
spreadsheet.  The  spreadsheet  contains  all  the  data  for  the  six  tracers, 
including  tritium,  carbon-14,  benzene,  p-xylene,  naphthalene,  and 
o-dichlorobenzene.  These  data  were  used  to  compare  measured  data  versus 
computed  data  at  various  locations  at  various  times. 
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3  Model  Development  with 
FEMWATER 


During  the  course  of  this  project,  the  MADE-2  modeling  system  was 
continuously  refined.  Two  numerical  grids  were  developed  and  several 
hydraulic  conductivity  fields  were  created.  Initially,  one  steady-state  flow 
field  was  used  to  represent  the  entire  test.  This  flow  field  represented  the 
average  water  table  condition,  which  is  roughly  approximated  by  the  water 
level  sampling  round  of  8  February  1991.  This  approximation  was 
improved  with  the  development  of  two  steady-state  flow  fields.  The  first 
flow  field  represents  the  low  water  table  conditions  during  the  first  half  of 
the  experiment.  The  second  flow  field  represents  the  high  water  table 
conditions  present  during  the  latter  half  of  the  experiment.  This  chapter 
discusses  FEMWATER  model  development  for  the  site,  beginning  with  a 
discussion  of  how  certain  features  of  the  software  system  were  modified 
during  this  work  to  make  the  GMS/FEMWATER  (Version  1.2)  software 
system  readily  applicable  to  sites  as  complex  as  the  CAFB  test  site. 


Summary  of  FEMWATER  Model 

FEMWATER  is  a  combination  of  two  public  domain  models, 
3DFEMWATER  (flow)  and  3DLEWASTE  (transport).  FEMWATER  may  be 
used  to  solve  flow  and  transport  of  chemicals  in  a  subsurface  using  the 
finite  element  numerical  technique.  Some  features  of  the  model  are  as 
follows: 

a.  It  treats  heterogeneous  and  anisotropic  media. 

b.  It  considers  both  distributed  and  point  source/sinks  that  are  spatially 
and  temporally  dependent. 

c.  It  accepts  a  prescribed  initial  condition  and  time-dependent  boundary 
conditions. 
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d.  It  automatically  determines  variable  boundary  conditions  of 
evaporation,  infiltration,  or  seepage  on  the  soil-air  interface. 

e.  It  provides  three  solution  options  for  estimating  the  nonlinear  matrix. 

/.  It  automatically  resets  time-step  size  when  boundary  conditions  or 
source/sinks  change  abruptly. 

g.  It  checks  the  mass  balance  computation  over  the  entire  region  for 
every  time-step. 


The  major  transport  processes  in  FEMWATER  are  advection, 
dispersion/diffusion,  adsorption,  first  order  decay,  and  source/sink.  For 
detailed  information  on  FEMWATER,  the  reader  is  referred  to  Lin  et  al. 
(1996). 


Modifications  to  the  GMS/FEMWATER  System 


There  are  a  number  of  features  in  the  GMS/FEMWATER  (Version  1.2) 
modeling  system  that  had  to  be  modified  during  this  work.  Some  changes 
were  made  to  the  FEMWATER  code  to  allow  this  model  application,  while 
other  features  will  require  further  modification.  A  discussion  of  how  the 
applied  modeling  approach  was  affected  or  altered  from  the  preferred 
course  follows. 


Maximum  allowable  material  properties 

The  MADE-1  and  MADE-2  experiments  were  designed  to  assess  the 
effect  of  heterogeneous  media  on  macrodispersion.  Much  work  has  been 
done  on  the  characterization  of  this  heterogeneous  media.  By  its  very 
nature,  the  porous  medium  at  the  site  is  highly  heterogeneous  and  this 
variability  causes  dispersion  of  any  tracer  moving  through  it.  Heterogeneity 
must  be  accounted  for  in  any  model  of  the  MADE  experiments. 

Since  Version  1.2  of  GMS/FEMWATER  allows  the  use  of  only  20 
material  properties,  this  means  only  20  hydraulic  conductivities  and  20  soil 
moisture  characteristic  curves  (which  define  porosity)  can  be  used.  With 
this  constraint,  achieving  a  highly  heterogeneous  model  was  not  possible 
without  some  modification  of  GMS/FEMWATER.  For  the  MADE 
experiment,  this  condition  had  to  be  removed.  All  necessary  arrays  in 
FEMWATER  (Version  1.2)  were  changed  so  that  the  number  of  hydraulic 
conductivities  that  could  be  specified  was  set  equal  to  the  number  of 
elements  used  in  the  model.  However,  Version  1.2  of  GMS  could  not  load  a 
large  number  of  material  types  without  an  extreme  degradation  in 
performance.  Once  GMS  (1.2)  reads  in  more  than  a  few  hundred  material 
properties,  GMS  begins  slowing  down  and  eventually  stops  reading  the 
materials  due  to  an  insufficient  memory  error.  To  avoid  the  performance 
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degradation,  for  graphic  demonstration  only,  a  geometry  file  having  only 
one  material  property  is  used  with  GMS  on  all  postprocessing  of  the  highly 
heterogeneous  model  simulations.  For  the  flow  and  transport  simulations 
with  FEMWATER  which  were  done  outside  of  GMS,  all  defined  materials 
and  hydraulic  conductivities  were  used. 


Maximum  allowable  boundary  conditions 

GMS/FEM WATER  (Version  1.2)  has  an  additional  hard-wired  limitation 
on  the  boundary  condition  specifications.  There  were  a  maximum  of  1,000 
allowable  boundary  conditions  for  rainfall-seepage,  Neuman  fluxes,  and 
Cauchy  fluxes.  The  model  developed  for  this  project  required  boundary 
condition  array  sizes  greater  than  1,000.  The  FEMWATER  code  was 
changed  so  that  up  to  3,000  boundary  conditions  could  be  specified. 


Conceptual  Groundwater  Flow  and  Transport 
Model  (FEMWATER) 


Hydrogeologic  conceptualizations  generally  must  consider  regional 
geology  well  beyond  the  extent  of  a  particular  study  site  (in  this  case,  the 
CAFB  test  site)  (Figure  2).  Regional  flow  trends  and  model  boundary 
conditions  must  be  defined.  Since  the  geological  setting  of  the  CAFB  site  is 
characterized  by  complex  fluvial  deposits,  it  is  important  to  understand  how 
these  deposits  are  distributed  and  how  they  affect  both  regional  and  local 
groundwater  flow  directions  and  rates. 

The  fluvial  deposits  consist  of  point  bar  deposits,  main  channel  deposits, 
channel  bottom  gravel  deposits,  and  overbank  deposits.  In  meandering 
systems  such  as  those  underlying  CAFB  (Figure  2),  complexities  arise  from 
channel  migration  and  avulsion  commonly  forming  meander  (oxbow) 
channels,  which  subsequently  fit  any  fine-grained  sediments.  The  MADE-2 
tracers  were  injected  into  finer-grained  (low  hydraulic  conductivity) 
overbank  deposits,  traveled  across  a  meander  channel  fill,  and  then 
reentered  overbank  deposits  that  are  mixed  with  sand  and  have  higher 
hydraulic  conductivity  (Figure  2). 

The  regional  groundwater  flow  direction  is  from  the  south  to  the  north, 
towards  the  Buttahatchee  River  (Figure  1).  The  northward  horizontal 
hydraulic  gradient  remains  unidirectional  throughout  the  tracer  test  area. 

The  vertical  hydraulic  gradients  vary  depending  upon  position  (Figure  4). 
Upward  gradients  exist  as  horizontal  water  movement  approaches  oxbow 
channels  because  of  the  increased  hydraulic  conductivity  (K)  of  these 
channels.  In  the  channels  themselves,  groundwater  flow  is  essentially 
horizontal,  with  few  to  no  vertical  gradients  existing.  Downward  vertical 
gradients  occur  down-gradient  of  the  channels  as  the  water  moves  into 
overbank  deposits  having  lower  hydraulic  conductivity.  During  transient 
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flow  events,  vertical  gradients  are  larger  than  during  low  water  table 
conditions. 

Flow  conditions  during  the  experiment  can  be  approximated  with  two 
steady-state  conditions.  During  the  first  170  to  180  days  of  the  experiment, 
no  recharge  affected  the  water  level  elevations,  as  shown  by  the  continuous 
records  in  Figure  5.  Subsequently,  numerous  water  level  responses  occur. 
Water  levels  in  the  shallow  fluvial  deposits  are  quite  variable  and  show 
rapid  response  to  the  rainfall.  The  continuous  records  indicate  that  water 
level  variations  of  up  to  4.7  m  have  occurred  during  the  440-day  experiment 
(Figure  5).  The  average  variation  during  this  period  from  all  16  wells  was 
3.9  m,  with  a  range  from  4.7  m  to  3.2  m.  When  the  eight  shallow  and  eight 
deep  wells  are  plotted  simultaneously,  it  is  clear  that  the  shallow  wells 
(Figure  6)  respond  with  the  same  timing  and  similar  water  level  responses 
as  the  deep  wells  (Figure  7).  This  indicates  that  direct  recharge  causes  the 
water  table  to  rise  in  the  aquifer.  Horizontal  hydraulic  gradient  reversals  do 
not  occur,  although  the  gradients  are  variable  because  spatial  variations  in 
porosity  and  hydraulic  conductivities  can  cause  different  water  level 
responses.  For  example,  the  horizontal  hydraulic  gradient  between  wells 
P53A  and  P61A  (Figure  3)  varies  from  0.0019  during  a  low  water  table 
condition  on  7  November  1990  to  0.0036  during  a  high  water  table 
condition  on  20  May  1991.  In  other  wells,  localized  horizontal  gradients 
can  increase  from  their  low-water  table  values  by  a  factor  of  three  to  four 
during  the  period  of  high  recharge  that  occurs  in  the  latter  half  of  the 
experiment. 

During  the  first  180  days  of  the  experiment,  the  groundwater  system  did 
not  receive  much  recharge  and  the  water  levels  showed  a  continuous 
decline.  This  period  of  slow  draining  of  the  aquifer  is  considered  a  near 
steady-state  condition.  Another  flow  condition  developed  thereafter  and 
continued  roughly  for  another  180-day  period,  followed  by  another  period 
of  declining  water  levels.  It  is  apparent  that  during  the  late  fall,  winter,  and 
into  early  spring,  the  groundwater  system  received  recharge.  During  the 
late  spring,  through  summer,  and  into  early  fall,  the  groundwater  system 
generally  drained.  The  tracer  test  was  affected  by  water  level  variations 
during  the  transient  conditions  for  two  reasons.  First,  the  increasing 
horizontal  hydraulic  gradients  during  the  transient  period  caused  faster 
tracer  movement.  Secondly,  the  additional  recharge  supplied  to  the  aquifer 
caused  significant  dilution  and  dispersion  of  the  leading  edge  of  the  tracers. 
This  aspect  of  the  conceptual  model  is  a  very  important  factor  to  consider 
when  the  tracer  test  concentration  data  are  analyzed. 


Areal  and  vertical  discretization  of  the  model 

The  plan  view  area  covered  by  this  model  is  21.0  x  103  m3  (5.2  acres) 
(Figure  4).  The  topographic  elevation  in  the  test  area  is  about  65  m  above 
mean  sea  level.  The  finite-element  nodes  were  stacked  13  high  with  each 
node  vertically  separated  by  1  m  (Figure  8).  The  nodal  layers  spanned  the 
elevations  from  65  m  down  to  53  m.  The  grid  is  variably  discretized  in  plan 
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Figure  5.  Continuous  well  hydrograph  of  P-53A,  which  is  located  at  coordinates  (-4.92,  -10.33)  or  11.5  m  southwest  of  the  injection  site 
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Figure  7.  Well  hydrographs  of  eight  continuously  recorded  “B-series”  (deep)  wells 
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Figure  8.  Orthographic  projection  of  the  three-dimensional  finite-element  grid 


view,  with  a  minimum  nodal  spacing  of  about  1  m  in  the  injection  area, 
enlarging  gradually  to  a  maximum  average  spacing  of  8  m.  The  total 
number  of  nodes  used  to  cover  this  area  is  10,621  and  the  total  number  of 
triangular  prismatic  (pentahedral)  elements  is  18,144,  The  total  volume  of 
the  model  is  275,275  m3.  Since  transport  through  the  saturated  zone  was  of 
primary  concern  in  this  study,  vertical  refinement  of  nodes  located  in  the 
vadose  zone  (  from  the  surface  downward  to  the  water  table)  was  not 
studied.  If  transient  flow  modeling  were  to  be  conducted,  finer  vertical 
discretization  would  be  required. 


Hydraulic  conductivity  distributions 

As  described  earlier,  the  GMS/FEMWATER  (Version  1.2)  program, 
under  review  here,  has  a  limit  of  20  material  properties.  For  an  application 
such  as  this,  20  hydraulic  conductivities  is  not  enough  considering  the 
degree  of  hydraulic  conductivity  variation  measured  with  the  borehole 
flowmeter.  Therefore,  for  this  application,  the  limit  of  20  material 
properties  was  removed  in  FEM WATER  (Version  1.2)  so  that  each  element 
could  have  its  own  hydraulic  conductivity.  This  precluded  the  use  of  GMS 
for  displaying  the  material  properties  because  there  are  18,144  hydraulic 
conductivities,  which  could  not  be  loaded  and  displayed  by  GMS  (1.2). 

The  hydraulic  conductivity  field  generated  for  this  model  used  the 
borehole  flowmeter  hydraulic  conductivities  available  in  the  MADE-2 
database.  A  total  of  2,018  individual  hydraulic  conductivities  (K)  were 
available  in  the  data  set.  For  each  borehole  location,  the  K  data  was 
vertically  averaged  over  1-m  intervals,  thereby  reducing  the  2,018  points  to 
405  points.  Vertical  averaging  was  done  because  the  vertical  element 
spacing  used  in  the  model  is  1  m.  The  range  of  measured  K’s  is  from  0.08 
to  1,249  m/day,  with  an  arithmetic  average  of  32.8  m/day. 

Using  the  vertically  averaged  hydraulic  conductivity  data,  three  main 
hydrogeologic  zones  consistent  with  the  fluvial  facies  model  may  be 
delineated  (Figure  3).  For  each  of  these  zones,  an  average  vertical  profile  of 
hydraulic  conductivity  was  determined  by  using  all  available  wells  in  the 
zone.  The  first  zone,  termed  the  near-field  overbank  deposits,  has  an 
overall  average  K  of  4.2  m/day  and  an  effective  porosity  of  0.31.  The 
second  zone,  termed  the  channel  deposits,  has  an  overall  average  K  of 
61  m/day  and  an  effective  porosity  of  0.11.  The  lower  effective  porosity 
seems  contradictory  considering  that  these  are  the  channel  deposits. 
Analysis  of  pumping  tests  indicates  that  three-dimensional  connections  of 
hydraulic  conductivity  exist,  with  relatively  narrow  high  conductivity 
channels  being  interconnected  in  three-dimensional  space.  These  channels 
allow  much  of  the  flow  through  the  aquifer,  with  an  estimate  that  over 
80  percent  of  the  aquifer  material  allows  less  than  10  percent  of  the  total 
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flow  through  the  aquifer.1  This  information,  combined  with  the  observation 
of  the  extreme  acceleration  of  the  leading  edge  of  the  plume  through  the 
channel  deposits,  indicates  that  these  deposits  have  an  overall  lower 
effective  porosity.  The  third  zone,  termed  the  far-field  overbank  deposits 
(intermixed  with  channel  sands),  has  an  overall  K  of  17  m/day  with  a  0.31 
effective  porosity.  Within  each  zone,  the  K  profile  changes  with  depth 
(Figure  8).  If  GMS  3-D  interpolation  is  completed  using  only  the  measured 
borehole  flowmeter  data,  the  interpolation  routines  have  to  extrapolate  out 
to  the  boundaries  of  the  model  because  flowmeter  data  do  not  exist  outside 
the  model  boundaries.  This  extrapolation  generally  created  hydraulic 
conductivities  which  are  heavily  weighted  towards  the  highest  or  lowest 
measured  data.  When  this  occurs,  the  hydraulic  conductivity  field  does  not 
fit  the  geologic  model  because  low  K  or  high  K  zones  are  erroneously 
(numerically)  created. 

To  avoid  this  problem,  randomly  located,  randomly  generated  K 
boreholes  were  placed  throughout  the  three  K  zones  so  that  the  GMS 
interpolation  routine  would  have  enough  information  throughout  the  model 
domain  to  avoid  extrapolation.  Thirty-five  additional  boreholes  were 
placed  in  zones  2  and  3,  while  seventeen  boreholes  were  placed  in  zone  1. 
The  vertical  K  profile  was  randomly  generated  using  the  average  K  vertical 
profile  for  each  of  the  three  zones.  The  randomness  within  each  borehole  at 
each  elevation  was  proportional  to  the  standard  deviation  in  measured  K’s. 
This  technique  creates  vertical  K  variations  in  the  boreholes  that  mimics  the 
measured  data.  With  the  additional  87  randomly  located  boreholes  added  to 
the  measured  data,  GMS  successfully  interpolated  the  3-D  K  field  to  the 
nodal  points  of  the  finite-element  domain.  Since  FEMWATER  (Version 
1.2)  requires  that  K  be  assigned  for  elements,  an  averaging  program  was 
written  that  took  the  six  nodal  K’s  for  each  element  and  arithmetically 
averaged  them.  When  viewed  with  a  3-D  graphics  program,  the  resulting  K 
field  fits  the  site  geologic  model  and  does  not  have  artificially  high  or  low 
K  zones. 

During  calibration,  no  changes  were  made  to  the  individual  hydraulic 
conductivities.  It  is  possible  that  the  model  calibration  could  be  improved 
by  manually  changing  the  K  field.  This  level  of  detail  was  not  attempted  in 
this  effort. 


Boundary  conditions 

There  are  essentially  five  sides  to  this  model,  including  a  top  and  bottom 
side,  and  three  lateral  sides.  Both  constant  head  and  flux  conditions  were 
used  in  the  steady-state  flow  models  (Figure  9),  whereas  interior  sources 
were  used  in  the  injection  area  to  introduce  the  concentration  for  the 
transport  models. 


5  Personal  communication  on  the  history  of  work  conducted  at  the  MADE  site,  1996,  Steve 
Young,  Environmental  Consulting  Engineers,  Inc.,  Knoxville,  TN. 
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Figure  9.  Orthographic  projection  showing  boundary  conditions  used  for  steady-state  flow 


General  description  of  flow  boundary  conditions 

Since  two  steady-state  flow  situations  were  simulated,  the  boundary 
conditions  used  represent  low  and  high  average  water  level  conditions 
during  the  experiment.  The  actual  water  level  data  sets  used  to  set  boundary 
conditions  were  measured  on  17  September  1990  (low  conditions)  and  8 
March  1991  (high  conditions).  Sensitivity  analysis  was  used  to  determine 
the  up-gradient  flux  boundary  conditions  needed  to  obtain  the  proper 
hydraulic  head  conditions  within  the  aquifer  for  each  time  simulated.  The 
down-gradient  boundary  conditions  were  set  to  constant  head  conditions, 
with  individual  boundary  conditions  set  for  each  time  simulated.  The 
constant  heads  set  along  the  down-gradient  boundaries  are  variable,  with  the 
highest  heads  on  the  east  side  of  the  model  and  the  lowest  on  the  west. 
Portions  of  the  lateral  boundaries  were  considered  no-flow  boundaries  and 
portions  were  considered  up-gradient  flux  conditions  (Figure  9).  A  uniform 
recharge  was  specified  through  the  top  of  the  model.  A  no-flow  condition 
was  specified  along  the  bottom. 

The  top  boundary  consisted  of  a  constant  flux  condition  to  represent 
recharge  to  the  aquifer.  The  constant  flux  was  set  equal  to  an  average 
recharge  rate  of  18.3  cm/year  (7.2  in. /year).  Since  no  rainfall  data  were 
available  in  the  MADE-2  database,  it  was  not  possible  to  determine  what 
percentage  of  rainfall  this  represents.  A  sensitivity  analysis  was  conducted 
to  investigate  the  importance  of  this  boundary  condition. 

A  uniform  no-flux  bottom  boundary  condition  was  used  at  53  m.  The 
transport  of  all  tracers  occurs  above  this  lower  boundary  throughout  the 
duration  of  the  experiment. 

The  lateral  boundaries  used  in  this  model  nearly  form  a  triangular  area 
(Figure  9).  The  apex  of  the  triangle  is  located  on  the  up-gradient  side  of  the 
model,  near  the  injection  area.  In  this  area,  horizontal  groundwater  flux  is 
allowed  to  enter  the  model,  to  simulate  the  regional  flow  of  water  into  the 
model  domain.  The  lateral  flux  used  along  this  boundary  is  either  constant 
from  top  to  bottom  (low  water  table)  or  is  vertically  variable  (high  water 
table).  This  boundary  condition  is  variable  because  during  transient 
groundwater  flow,  larger  vertical  gradients  exist  along  this  boundary  and  a 
vertically  variable  flux  is  needed  to  generate  these  gradients.  For  the  high 
water  table  conditions,  the  up-gradient  fluxes  are  greatest  near  the  water 
table  and  decrease  downward.  The  average  of  these  fluxes  is  equal  to  the 
measured  hydraulic  gradient  across  the  boundary  multiplied  by  the  average 
hydraulic  conductivity  of  the  up-gradient  area.  In  other  words,  the  fluxes 
are  physically  realistic. 

The  down-gradient  boundary  is  set  to  a  constant  head.  Since  this 
boundary  is  at  an  angle  with  respect  to  the  regional  flow  system,  variable 
constant  heads  are  used  along  its  length  (Figure  9).  Hydrostatic  heads  are 
set  for  each  vertical  column  of  nodes  and  these  head  conditions  are  based  on 
the  measured  down-gradient  water  levels  each  time  period  simulated. 
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General  description  of  transport  boundary  conditions 

The  transport  boundary  conditions  were  set  to  represent  injection  of  the 
tracers.  The  maximum  measured  concentration  from  each  sampling  period 
was  determined  and  plotted  with  respect  to  time.  These  data  indicated  that 
the  concentrations  generally  decreased  over  time,  as  would  be  expected. 
From  these  data,  transient  concentrations  were  set  at  the  injection  well. 
Numerous  simulations  were  performed  using  various  boundary  conditions 
for  the  injection  wells,  and  this  decreasing  concentration  condition  was 
found  to  be  the  best  representation  of  the  experiment.  The  use  of  constant 
concentrations  at  the  injection  well  was  shown  to  be  inappropriate  through 
numerous  simulations.  No  other  transport  boundary  conditions  were 
required  because  the  plume  stayed  within  the  model  boundaries  during  the 
entire  simulation. 

Interior  sources  were  used  to  represent  the  injection  of  tracer  into  the 
model.  Four  vertically  stacked  nodes  at  coordinates  (0.5,  0.3  m)  were  used 
to  represent  the  injection  interval.  These  nodes  span  the  elevation  of  61  m 
to  57  m.  The  actual  screened-section  openings  in  the  injection  wells  span 
the  range  of  58.1  to  57.5  m.  Due  to  increased  fluid  pressure  in  the  injection 
well  area,  the  injected  concentrations  were  initialized  slightly  above  and 
below  the  screened-section  interval. 
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4  FEMWATER  Application 
and  Results 


The  modeling  approach  adopted  here  simulates  the  experiment  with  a 
combination  of  two  steady-state  flow  fields  and  transient  transport.  The 
first  flow  field  represents  early-time,  low-water-table  conditions  from  day  0 
(starting  day)  to  day  180.  The  second  flow  field  represents  high-water-table 
conditions  during  the  latter  portion  of  the  test  (day  180  to  450).  During  this 
GMS/FEMWATER  application,  it  was  determined  that  transient  flow  and 
transient  transport  would  require  too  much  time  to  successfully  complete. 
The  computational  demands  of  that  approach  would  substantially  limit  the 
number  of  modeling  runs  that  could  be  completed. 

During  the  course  of  this  study,  several  iterations  on  the  model  grid 
discretization  and  on  representation  of  the  three-dimensional  hydraulic 
conductivity  field  were  made.  Once  these  were  defined,  numerous  model 
runs  were  made  to  determine  the  model  sensitivity  to  the  applied  boundary 
conditions.  The  best  boundary  conditions  were  selected  for  both  the  low 
and  high  water  table  conditions,  and  steady-state  solutions  were  achieved. 
These  solutions  are  termed  calibrated  conditions,  but  these  really  represent 
a  reasonable  comparison  between  the  calculated  and  measured  water  levels 
and  hydraulic  gradients  across  the  entire  site. 

To  ensure  that  FEMWATER  source  code  changes  made  during  this 
project  were  valid,  code  accuracy  had  to  be  verified  after  the  code  changes 
were  completed.  A  simple  two-dimensional  flow  and  transport  example 
was  created  so  that  the  FEMWATER  flow  and  transport  results  could  be 
compared  to  the  analytical  solution  for  this  problem.  This  example  is  a 
simple  tracer  experiment  in  which  a  constant  concentration  is  injected  along 
a  boundary  into  a  uniform  flow  field.  The  growth  of  the  example  plume  is 
dependent  upon  the  ratio  between  the  velocity  and  dispersion  terms  used. 
Although  a  detailed  numerical  analysis  was  not  conducted  to  compare  the 
results  of  the  example,  graphical  comparisons  of  results  appeared  to  be 
reasonable. 
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Use  of  Gauge  Plots  and  Methods  of  Model 
Calibration 


One  of  the  goals  of  site-specific  modeling  is  to  compare  measured 
parameters  to  simulated  parameters.  To  accomplish  this,  GMS  allows 
comparison  of  measured  values  at  one  well  (gauge)  to  the  model-calculated 
values  at  the  same  location.  Up  to  five  comparisons  can  be  made  at  any  one 
location  so  that  water  level,  concentration,  velocity,  etc.  can  be  compared. 

It  is  also  expected  in  GMS  that  measurements  will  have  been  made  at  more 
than  one  time  (i.e.,  data  are  transient).  This  presents  two  difficulties  that 
have  to  be  overcome.  First,  the  comparisons  can  be  made  at  only  one 
location  at  a  time.  This  type  of  comparison  works  best  for  a  tracer 
breakthrough  curve,  for  example,  but  it  is  not  efficient  for  comparing  a 
complete  water  level  round  to  the  calculated  values.  Secondly,  it  is  not 
possible  to  compare  single  readings  at  any  location,  such  as  one  might 
generate  during  a  water  table  elevation  sampling  round.  Therefore,  the  use 
of  gauge  plots  makes  it  very  difficult  to  determine  how  well  a  flow  model  is 
calibrated  at  a  particular  time.  A  common  calibration  practice  is  to  compare 
a  set  of  water  levels  measured  at  numerous  wells  to  model-calculated  water 
levels  (Figure  10),  for  any  given  time.  This  type  of  information  is  valuable 
in  calibrating  a  model  because  well-by-well  comparisons  can  be  made  and 
inaccuracies  in  the  model  can  be  revealed  quickly.  Having  a  software  for 
calibration  purposes  in  numerical  groundwater  flow  models  is  very  useful 
and  allows  a  rapid  analysis  of  the  well-  by-well  water  level  comparisons. 
Normally  a  file  consisting  of  well  coordinates  and  measured  water  levels  at 
various  times  is  created  from  the  measured  data.  This  file  is  used  by  the 
model  after  convergence  is  reached  so  that  interpolations  of  the  model- 
calculated  water  levels  can  be  compared  to  the  measured  data.  Statistical 
summaries  of  the  comparison  can  be  created  and  used  to  assess  or  document 
how  well  a  model  is  calibrated. 

GMS/FEMWATER  (Version  1.2)  did  not  allow  this  type  of  comparison 
to  be  made.  Therefore,  it  is  difficult  to  assess  how  well  a  model  is 
calibrated.  The  gauge  plot  methodology  in  GMS  (Version  1.2)  requires  that 
each  well  be  compared  individually  to  determine  how  well  model-computed 
water  levels  compare  to  measured  values.  In  applied  modeling  this  is 
inefficient  and  a  change  needs  to  be  made  so  that  a  simultaneous 
well-by-well  comparison  can  be  made.  For  this  reason,  in  this  application, 
computed  versus  measured  water  levels  are  compared  graphically. 


Steady-State  Flow  Calibration 

As  with  any  model,  the  solutions  obtained  for  the  hydraulic  head  field 
can  be  achieved  with  different  sets  of  input  parameters  because  the 
solutions  are  not  unique.  The  intent  of  this  modeling  was  not  to  make  the 
most  highly  calibrated  model  possible,  but  rather  to  create  a  model  that  used 
as  much  of  the  available  site  data.  In  light  of  this  philosophy,  the 
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Figure  10.  Example  plot  comparing  model  calculated  to  measured  hydraulic  heads  (data  not  from  MADE-2  model) 


conductivities  were  not  altered  on  an  element-by-element  basis  to  improve 
the  calibration,  nor  was  variable  recharge  applied  in  any  areas.  The  final 
hydraulic  conductivity  field  used  is  based  only  on  measured  data. 

The  primary  goal  in  developing  the  low-  and  high-water-table  models 
was  to  create  two  hydraulic  gradient  fields  that  could  be  used  to  advect  the 
tracers.  These  gradient  fields  have  different  properties,  with  the  average 
gradient  being  higher  during  the  high-water-table  condition.  During 
low-water-table  conditions,  vertical  gradients  are  minimal  and  the 
horizontal  hydraulic  gradient  through  the  injection  area  is  less  than  during 
high-water-table  conditions.  The  increased  hydraulic  gradients  that  exist 
during  the  latter  half  of  the  experiment  cause  the  tracers  to  move  through 
the  channel  faster.  In  addition,  the  recharge  that  occurs  during  the  second 
half  of  the  experiment  dilutes  and  disperses  the  tracers  that  exist  in  the 
channel  deposits. 

The  average  hydraulic  head  difference  across  the  site  during  the 
low-water-table  condition  is  0.72  m  and  is  0.95  m  during  high-water-table 
conditions.  The  total  head  steady-state  solutions  for  the  low-water-table 
condition  (Figure  11)  and  the  high-water-table  condition  (Figure  12)  are 
shown  along  with  the  measured  water  levels.  The  water  levels  plotted  at  the 
well  locations  represent  either  a  single  water  level  reading  taken  at  a  well  or 
vertically  averaged  values,  if  a  nested  well  exists  at  that  location.  As  can  be 
seen  in  these  plots,  differences  exist  between  calculated  and  measured  water 
levels  (total  heads)  throughout  the  model  domain.  Once  again,  no  attempt 
was  made  to  reduce  the  difference  between  computed  and  simulated  water 
levels.  The  primary  goal  was  to  obtain  a  reasonable  approximation  to  the 
water  table  in  general  and  to  make  sure  the  hydraulic  gradients  were  correct 
in  magnitude  and  direction. 

The  velocity  fields  corresponding  to  these  solutions  (Figures  13  and  14) 
generally  show  unidirectional  flow,  from  the  south  to  the  north.  The 
velocity  profiles  vary  more  for  the  high-water-table  condition  because  the 
up-gradient  flux  boundary  condition  was  vertically  variable,  whereas  it  is 
uniform  during  low-water-table  conditions.  The  average  along-gradient 
velocity  in  the  high-water-table  condition  is  2.2  times  greater  than  during 
the  low-water-table  condition.  The  vertical  variation  in  velocity  is  also 
greater. 

These  flow  solutions  were  used  as  input  to  the  tritium  transport  model. 
The  low- water-table  solution  was  used  for  the  first  180  days  of  the 
experiment.  The  concentration  field  solution  for  day  180  was  used  as  an 
initial  condition  for  transport  during  the  high-water-table  condition. 
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Figure  13.  Orthographic  projection  of  steady-state  Darcy  velocity  for  the  low-water-table  condition 
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Transport  Analysis 


Simulation  of  the  tritium  plume  was  completed  for  the  entire  experiment. 
Numerous  runs  were  made  to  test  the  sensitivity  of  the  model  to  various 
parameters  including  the  dispersion  terms,  the  particle  tracking  parameters 
of  FEM WATER,  time-step  selection,  and  effective  porosity. 

Two  terms  representing  longitudinal  and  transverse  dispersivity  are 
available  in  FEMWATER.  Simulations  using  longitudinal  dispersivity 
ranging  from  10  m  to  0.5  m  were  made.  Lateral  dispersivities  ranging  from 
5  m  to  0.01  m  were  tested.  In  all  cases,  if  any  dispersion  terms  are  used,  the 
computed  concentrations  will  move  up-gradient  of  the  injection  area,  even 
though  no  up-gradient  velocity  terms  exist.  The  larger  the  dispersion  term 
used,  the  more  that  up-gradient  plume  spreading  occurs.  This  up-gradient 
spreading  occurs  even  if  a  very  small  time-step  is  used.  Short-term  models 
were  run  with  time-steps  of  0.001  day,  and  up-gradient  spreading  occurred. 
This  up-gradient  spreading  occurred  even  when  the  particle  tracking 
parameters  ( nx ,  ny,  and  nz)  in  GMS/FEMWATER  were  increased  from  2,2,2 
to  4,4,4.  These  terms  are  used  to  split  elements  in  each  direction  to  improve 
the  particle  tracking  accuracy.  Although  difficult  to  quantify,  it  seemed  that 
by  increasing  these  terms,  less  up-gradient  spreading  might  occur. 

However,  model  run  time  increased  because  the  time  for  particle  tracking 
increased. 

With  a  highly  heterogeneous  hydraulic  conductivity  field  in  this  model, 
inherent  numerical  dispersion  occurs  in  these  simulations.  Tortuous  flow 
paths  can  occur  because  different  hydraulic  conductivities  are  found  from 
element  to  element.  Computed  versus  measured  concentrations  are  shown 
in  series  of  four  representing  the  sampling  rounds  of  27,  132,  224,  and  328 
days.  These  rounds  were  conducted  about  100  days  apart.  With  an 
injection  concentration  of  55,610  pCi/ml,  the  tritium  results  are  shown  for 
50  pCi/ml  (10 ‘3  C/Co)  and  5  pCi/ml  (10"4C/Co).  The  computed  plumes  are 
shown  as  blue  isosurfaces,  whereas  the  measured  data  are  shown  as  a  solid 
circle,  square,  or  triangle.  These  point  measurements  represent  every 
location  where  a  measurement  was  taken  and  the  value  was  greater  than  the 
contoured  value.  For  example,  for  a  50-pCi/ml  plot,  every  solid  point 
represents  a  reading  of  at  least  50  pCi/ml.  The  sampling  points  that 
measured  less  than  50  pCi/ml  are  not  shown.  The  50-pCi/ml  plumes  for 
days  27,  132,  224,  and  328  are  shown  in  plan  view  (Figures  15-18)  and 
vertical  view  (Figures  19a-19d).  From  an  examination  of  these  plots, 
several  features  of  the  transport  simulation  can  be  determined.  First,  the 
areal  extent  of  the  27th  day  plume  is  underpredicted  because  the  effect  of 
the  two-day  injection  is  not  included  in  the  simulation.  The  early  time 
plume  is  larger  than  would  be  expected,  given  the  flow  field  because  the 
injection  forced  the  tracers  to  move  in  directions  other  than  just 
down-gradient.  Secondly,  the  50-pCi/ml  plume  is  simulated  quite  well 
through  the  224-day  period  (Figures  15-17,  19a, b,c).  This  sampling  round 
is  about  40  days  into  the  transient  flow  portion  of  the  test.  This  implies  that 
the  first  portion  of  the  test,  the  low-water-table  condition,  is  simulated  quite 
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Figure  16.  Plan  view  of  50-pCi/ml  tritium  plume  at  132  days 
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Figure  19.  Vertical  view  of  50-pCi/ml  tritium  plume  at:  a)  27  days,  b)  132  days,  c)  224  days,  d)  328  days 
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well.  By  the  time  the  328-day  sampling  round  was  conducted,  the 
50-pCi/ml  plume  was  predicted  to  be  about  100  m  further  down-gradient 
than  it  was  measured  (Figure  18).  This  implies  that  enhanced  dilution  from 
the  recharge  during  the  transient  portion  of  the  test  diluted  the  leading  edge 
of  the  plume  such  that  no  readings  over  50  pCi/ml  were  measured.  This  is 
supported  by  the  fact  that  the  50-pCi/ml  plume  was  measured  further 
down-gradient  at  day  224  than  it  was  at  day  328.  Dilution  was  an  important 
process  operating  during  the  last  half  of  this  experiment. 

The  5-pCi/ml  plumes  for  days  27,  132,  224,  and  328  are  shown  in  plan 
view  (Figures  20  -  23).  The  down-gradient  plan  view  extent  of  this  plume 
is  always  underpredicted.  The  early  time  data  show  that  one  effect  of 
injection  was  to  move  the  5-pCi/ml  plume  about  25  m  further 
down-gradient  than  would  be  expected  due  to  advection.  If  25  m  were 
added  to  each  calculated  plume  extent,  the  simulations  for  days  27  and  132 
would  be  very  accurate.  The  plume  longitudinal  extents  for  days  224  and 
328  are  also  underpredicted,  but  this  time  by  more  than  25  m.  The  lateral 
extents  of  these  plumes,  however,  are  very  accurate.  Lateral  dispersion  of 
this  plume  is  simulated  quite  well.  Without  additional  recharge  during  the 
transient  phases  of  this  test,  however,  it  is  not  possible  to  get  the  5-pCi/ml 
plume  to  accelerate  all  the  way  to  the  down-gradient  edge  of  the  model 
while  keeping  the  50-pCi/ml  plume  limited  to  a  down-gradient  extent  of 
125  m. 

To  add  another  piece  of  evidence  to  the  analysis  that  indicates  the 
pronounced  effect  of  dilution  after  day  180,  GMS  gauge  plots  can  be  used. 
Gauge  plots  are  used  to  show  measured  versus  computed  values  at  any 
location  in  the  model  domain.  Eight  locations  were  randomly  selected  from 
the  concentration  database  (Figure  24).  From  two  of  these  near-field 
locations,  tritium  break-through  curves  are  displayed  (Figures  25  -  26).  At 
location  1018,  which  is  located  at  coordinates  -0.1,13.5,58.85,  the  model 
predicts  the  timing  and  magnitude  of  the  tritium  peak  concentration  quite 
well  (Figure  25).  Between  the  sampling  round  of  days  132  and  224,  the 
observed  tritium  declined  from  1,987  to  1,940  pCi/ml.  During  this  time, 
this  sampling  point  was  located  between  2.87  and  5.78  m  beneath  the  top  of 
the  water  table.  For  the  remainder  of  the  sampling  events,  this  location 
showed  a  typical  tailing-edge  of  the  concentration  breakthrough  curve.  In 
other  words,  the  additional  recharge  did  not  cause  dilution  at  this  point, 
which  was  at  least  2.87  m  below  the  top  of  the  water  table. 

At  location  1644,  which  is  located  at  coordinates  -4.8,25.9,60.23,  the 
measured  tritium  breakthrough  curve  was  quite  different.  Between  the 
sampling  rounds  of  days  132  and  224,  the  concentration  dropped  from  73  to 
2  pCi/ml.  The  modeled  tritium  concentration  also  shows  the  same  effect. 
The  modeled  effect  occurred  because  at  day  180,  the  water  table 
instantaneously  rose  from  the  low-  to  the  high-water-table  condition.  This 
instantaneously  caused  dilution,  as  can  be  seen  in  the  breakthrough  curve. 
With  this  sampling  point  being  only  1.13  m  below  the  low  water  table 
position,  it  is  more  likely  that  the  effect  of  dilution  can  be  observed. 
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Figure  20.  Plan  view  of  5-pCi/ml  tritium  plume  at  27  days 


Figure  21 .  Plan  view  of  5-pCi/ml  tritium  plume  at  1 32  days 
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Figure  22.  Plan  view  of  5-pCi/ml  tritium  plume  at  224  days 
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Figure  23.  Plan  view  of  5-pCi/ml  tritium  plume  at  328  days 
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Observed  versus  computed  Tritium 


Figure  25.  Concentration  curve  at  location  1018  showing  measured-versus-computed  tritium  (pCi/ml) 
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Observed  versus  computed  Tritium 


Figure  26.  Concentration  curve  at  location  1644  showing  measured-versus-computed  tritium  (pCi/ml) 
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Summary  of  Transport  Results 


Although  two  steady-state  flow  fields  were  used  to  generate  the  tritium 
plume,  it  does  not  appear  possible  to  accurately  simulate  this  experiment 
with  steady-state  flow  fields.  After  extensive  model  testing  and  analysis  of 
the  concentration  data,  the  dilution  effects  caused  by  transient  recharge 
during  the  last  half  of  the  experiment  make  it  virtually  impossible  to  match 
the  low-concentration,  leading  edge  of  the  tracer  plumes.  For  example,  the 
tritium  plume  has  a  peculiar  characteristic  in  that  much  of  the  high 
concentration  portions  of  the  plume  stay  in  the  injection  area  during  most  of 
the  experiment.  During  the  entire  duration  of  the  experiment,  the  50-pCi/ml 
isopleth  (10'3  C/C  )  moved  a  maximum  distance  of  roughly  125  m 
down-gradient  of  the  injection  area.  The  total  transport  distance  is  over 
300  m  for  far-field  concentrations  of  5  pCi/ml  (10'4  C/CQ) .  This  disparity 
between  the  50-  and  5-pCi/ml  isopleths  makes  this  experiment  very  difficult 
to  simulate.  Transient  recharge  fluxes  and  hydraulic  gradients  are  necessary 
to  cause  dilution  and  dispersion  at  the  leading  edge  of  the  tritium  plume. 

The  effects  of  these  two  features  cannot  be  captured  within  steady-state 
flow  models. 


Sensitivity  Analysis 

The  flow  model  was  tested  for  sensitivity  to  the  up-gradient  horizontal 
flux  boundary  condition  and  the  recharge  rate.  The  transport  model  was 
tested  for  sensitivity  to  the  dispersion  terms,  effective  porosity,  and 
time-step  selection.  The  results  of  the  transport  sensitivity  have  already 
been  discussed  in  the  previous  section. 

Generally,  sensitivity  analysis  is  performed  on  a  model  after  it  is 
calibrated.  Changes  are  made  to  input  parameters  and  model  response  is 
evaluated.  In  this  case,  because  it  is  difficult  to  quantify  calibration  with 
Version  1.2  of  the  GMS/FEMWATER  software  system,  only  general  model 
behavior  can  be  reported.  The  new  version  of  GMS  (Version  2.0)  may  have 
improved  numerical  tools  for  model  calibration.  Sensitivity  analysis  was 
completed  on  the  important  input  parameters. 

The  up-gradient  horizontal  flux  boundary  condition,  along  with  direct 
recharge  through  the  top  of  the  model,  controls  water  levels  throughout  the 
model.  The  maximum  water  levels  recorded  during  the  experiment  always 
occur  along  the  up-gradient  boundary.  A  series  of  model  runs  were 
conducted  to  show  the  relationship  between  the  up-gradient  flux  boundary 
condition  and  the  maximum  water  levels  in  the  model  (Figure  27).  For 
these  runs,  recharge  was  held  constant  at  18.3  cm  (7.2  in.)  per  year.  It  can 
be  seen  that  a  nearly  linear  relationship  exists  between  the  up-gradient 
boundary  flux  and  the  calculated  heads  in  the  domain.  For  every  1-m  rise  in 
the  water  table,  the  up-gradient  flux  has  to  be  increased  by  a  factor  of  4. 

This  relationship  is  linear  throughout  the  entire  range  of  up-gradient  fluxes 
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Figure  27.  Maximum  total  head  versus  calibrated  up-gradient  flux  multiplier 


measured.  The  maximum  horizontal  velocity  corresponding  to  these  cases 
is  not  quite  linear,  however  (Figure  28).  If  the  up-gradient  flux  is 
multiplied  by  a  factor  greater  than  1,  the  maximum  velocity  also  increases 
by  that  factor.  If  the  multiplicative  factor  is  less  than  1,  the  velocity 
diminishes  less  rapidly.  In  this  case,  the  up-gradient  boundary  condition 
heads  are  nearing  the  down-gradient  heads,  and  the  relationship  begins 
breaking  down  because  up-gradient  flow  begins  occurring. 

Similarly,  a  series  of  model  runs  were  conducted  for  various  recharge 
rates,  holding  the  up-gradient  horizontal  flux  constant.  Results  from  these 
tests  indicate  that  the  overall  mass  balance  of  this  model  is  satisfied  by  the 
natural  balance  between  the  up-gradient  flux  entering  the  model  and  the 
down-gradient  boundary  flux  leaving  the  model.  No  matter  what  recharge 
flux  was  used,  the  recharge  contribution  was  only  a  small  fraction  of  the 
total  water  movement  through  the  model.  If  the  up-gradient  flux  was  not 
used,  the  entire  recharge  rate  was  correctly  calculated  by  FEMWATER.  If 
up-gradient  flux  was  used,  the  calculated  recharge  rate  did  not  equal  what 
was  used  in  the  boundary  conditions  —  it  was  always  less  than  expected. 
Tests  were  even  conducted  in  which  linear  soil  moisture  characteristic 
curves  were  used.  In  this  case,  slightly  more  recharge  entered  the  model 
but  it  still  did  not  equal  the  rate  specified  on  the  boundary  conditions. 

There  are  three  possible  explanations  for  this  behavior.  The  up-gradient 
boundary  condition  is  being  specified  along  a  triangular-shaped  model 
boundary  that  happens  to  be  located  in  a  zone  of  converging  flow.  Even 
with  analyzing  the  available  water  level  data  from  this  area,  it  is  difficult  to 
tell  how  much  of  the  up-gradient  boundary  should  be  set  as  a  flux  condition. 
The  areal  coverage  currently  specified  by  the  up-gradient  boundary  flux 
conditions  might  be  larger  than  required.  This  implies  that  more  water 
might  be  horizontally  moving  into  the  model  than  is  actually  occurring  in 
this  area.  The  net  effect  is  that  this  additional  water  overrides  the  effect  of 
recharge.  Secondly,  because  recharge  is  entering  through  the  unsaturated 
zone,  the  total  flux  is  being  computed  with  a  relative  hydraulic  conductivity 
term  that  is  less  than  one,  because  the  top  of  the  model  is  always 
unsaturated  in  the  steady-state  solutions.  It  is  possible  that  a  transient  flow 
model  would  produce  a  more  accurate  model  containing  recharge. 


Conclusions  on  FEMWATER  Application 

A  careful  examination  of  the  MADE-2  hydrogeologic  databases  is 
necessary  to  fully  realize  the  difficulty  in  numerically  simulating  this 
experiment,  especially  during  the  last  half  of  the  experiment.  Movement  of 
the  tracers  was  strongly  influenced  by  the  heterogeneous  porous  media  in 
which  the  experiment  took  place  and  by  the  transient  flow  system  during 
the  latter  half  of  the  experiment.  Simulated  results  of  the  first  180  days  of 
the  experiment  show  that  tracer  concentrations  can  be  reproduced  quite  well 
with  FEMWATER.  The  concentration  distributions  of  the  latter  half  of  the 
experiment  are  difficult  to  reproduce.  The  transient  rainfall  recharge  that 
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Figure  28.  Maximum  horizontal  Darcy  velocity  versus  calibrated  up-gradient  flux  multiplier 


occurred  during  this  time  had  two  net  effects  which  must  be  accounted  for 
to  obtain  accurate  numerical  results.  First,  the  hydraulic  gradients 
increased,  at  times  by  as  much  as  a  factor  of  4.  Secondly,  the  additional 
infiltrated  water  diluted  the  plume,  causing  reduced  and  dispersed 
concentrations  at  the  leading  edge  of  the  plume.  With  a  steady-state  flow 
field  approach,  these  effects  could  not  be  accurately  reproduced. 


Suggestions  for  Additional  Work 

Additional  modeling  could  attempt  to  perform  transient  flow  and 
transport  modeling  of  the  latter  half  of  the  experiment.  The  conditions  that 
will  be  required  for  success  in  this  endeavor  include  obtaining  a  successful 
dynamic  balance  between  three  transient  data  sets.  Transient  up-gradient 
fluxes,  down-gradient  head  profiles,  and  recharge  fluxes  must  be  balanced. 
The  transient  recharge  fluxes  will  have  to  be  derived  from  the  continuous 
water  level  records  measured  during  the  experiment.  If  a  rainfall  record  is 
available  for  the  experiment  (it  is  not  in  the  database),  the  rainfall-seepage 
boundary  condition  could  be  used  in  FEMWATER.  In  any  case,  it  is 
expected  that  transient  flow  and  transport  analysis  of  this  portion  of  the 
experiment  will  be  computationally  very  demanding.  The  numerical  grid 
will  require  additional  vertical  refinement  between  elevations  65  and  62  m, 
with  cells  placed  every  0.2  m.  The  new  grid  necessary  for  transient  flow 
would  then  have  at  least  25  vertical  layers,  with  greater  than  20,000  nodes 
and  nearly  40,000  elements. 
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5  MODFLOW  Modeling 
Approach 


The  remainder  of  this  report  describes  the  application  of  GMS/ 
MODFLOW-MT3D  to  the  CAFB  test  site.  In  MODFLOW/MT3D,  the 
modeling  domain  was  discretized  to  a  structured  (rectangular)  grid  system. 
The  distribution  of  hydrogeologic  data  input  for  the  grid  system  is  different 
than  the  mesh  system  used  by  FEMWATER.  For  a  detailed  description  of 
FEMWATER  and  MODFLOW/MT3D  data  requirements,  the  reader  is 
referred  to  the  GMS  user’s  manual  (Brigham  Young  University  1995).  The 
MADE  site  hydrogeologic  data  were  interpolated  again  to  fit  data 
requirements  for  MODFLOW/MT3D.  This  section  describes  the 
development  of  a  conceptual  model  for  application  of  MODFLOW/MT3D 
using  MADE-2  data  sets.  Descriptions  of  flow  and  transport  simulations 
have  been  separated  into  two  chapters,  Chapter  5  for  flow  and  Chapter  6  for 
transport. 


Summary  of  MODFLOW  Model 

MODFLOW  is  a  widely  used  three-dimensional,  finite-difference 
groundwater  (saturated  zone)  flow  model.  It  has  a  modular  structure  that 
allows  it  to  be  modified  for  a  particular  application.  Many  new  capabilities 
have  been  added  to  the  original  model.  The  GMS/MODFLOW  version  used 
in  this  demonstration  was  the  latest  one  available  from  the  U.S.  Geological 
Survey. 

The  model  simulates  steady  and  unsteady  flow  in  aquifer  layers  which 
can  be  confined,  unconfined,  or  a  combination  of  both.  Flow  from  external 
stresses,  such  as  flow  to  wells,  areal  recharge,  evapotranspiration,  flow  to 
drains,  and  flow  through  riverbeds,  can  be  simulated.  Hydraulic 
conductivities  or  transmissivities  for  any  layer  may  differ  spatially 
(heterogeneous)  and  be  anisotropic.  The  storage  coefficient  may  also  change 
spatially.  Specified  head  and  specified  flux  boundary  conditions  can  be 
imposed.  Head-dependent  flux  across  the  model’s  outer  boundary  can  be 
specified.  This  option  allows  water  to  be  supplied  to  a  boundary  block  in  the 
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modeled  area  at  a  rate  proportional  to  the  current  head  difference  between  a 
“source”  of  water  outside  the  modeled  area  and  the  boundary  block. 

The  groundwater  flow  equation  is  solved  using  approximated 
block-centered  finite  differences.  Several  solutions  are  provided  to  solve  the 
associated  matrix  problem;  the  user  can  choose  the  best  solution  for  the 
particular  problem.  Mass  balances  are  computed  for  each  time-step  and  as  a 
cumulative  volume  from  each  source  and  type  of  discharge. 

In  order  to  use  MODFLOW,  initial  conditions,  hydraulic  properties,  and 
stresses  must  be  specified  for  every  model  cell  in  the  finite-difference  grid. 

Primary  output  is  head,  which  can  be  written  to  the  listing  file  or  into  a 
separate  file.  Other  output  includes  the  complete  listing  of  all  input  data, 
drawdown,  and  budget  data.  Budget  data  are  printed  as  a  summary  in  the 
listing  file,  and  detailed  budget  data  for  all  model  cells  can  be  written  into  a 
separate  file.  The  hydraulic  head  distribution  can  be  visualized  by  GMS. 

For  detailed  information  on  MODFLOW,  the  reader  is  referred  to 
McDonald  and  Harbaugh  (1988). 


Model  grid 

The  final  grid  developed  and  used  with  GMS  for  the  site  study  consisted 
of  25,025  blocks  comprised  from  11  layers,  65  rows,  and  35  columns  as 
shown  in  Figure  29.  A  uniform  grid  with  block  dimensions  of  5  x  5  x  1  m 
was  used.  The  grid  spacing  used  is  relatively  coarse  for  the  degree  of 
heterogeneity  at  the  MADE  site,  but  was  chosen  because  of  computational 
speed  and  data  storage  requirements.  The  elevation  of  the  top  of  the  first 
layer  was  set  to  62.5  m  with  the  thickness  of  each  underlying  layer  being 
1  m  and  a  bottom  elevation  of  the  lowermost  layer  of  51.5  m. 


Hydraulic  conductivity 

The  available  hydraulic  conductivity  data  from  the  MADE  experiments 
include  2,187  measurements  (borehole  flowmeter  data)  for  67  wells 
(Figure  30).  Hydraulic  conductivity  varies  from  10"4  cm/s  to  10"2  cm/s.  The 
variance  of  the  natural  logarithm  of  hydraulic  conductivity  for  the  aquifer  is 
4.5  (Rehfelt,  Boggs,  and  Gelhar  1992;  Boggs  et  al.  1990).  This  high 
variance  indicates  that  the  site  is  one  of  the  most  heterogeneous  sites 
extensively  investigated  and  is  an  order  of  magnitude  larger  than  that  for 
the  Borden  and  Cape  Cod  test  sites,  both  of  which  have  measurement  scales 
comparable  to  the  MADE  site. 

The  hydraulic  conductivity  data  were  interpolated  to  estimate  values  at 
each  cell  of  the  grid.  GMS  has  several  data  interpolation  options  including 
inverse  distance  weighting,  nearest  neighbor,  and  kriging.  These  options 
were  applied  to  the  hydraulic  conductivity  data  and  the  results  were 
compared.  To  apply  field  hydrogeologic  conditions  into  the  numerical  grid, 
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the  2-D  and  3-D  data  interpolation  routines  available  under  GMS  1.2  and  a 
commercial  interpolation  package  were  used.  The  interpolated  hydraulic 
conductivity  numbers  were  incorporated  into  the  numerical  grid  using  the 
GMS  scalar  data  file  format. 


Piezometric  head 

Piezometric  head  data  were  measured  at  approximately  2-week  intervals 
during  MADE-2.  The  accuracy  of  the  measurements  is  approximately  ±0.6 
cm.  Data  were  available  for  manual  measurements  and  some  continuous 
monitoring  wells.  This  study  only  used  the  manually  measured  data 
because  more  wells  were  measured  manually  than  continuously.  Table  2 
below  indicates  the  dates  for  which  manual  piezometric  measurements  were 
taken  and  Figure  31  shows  the  location  of  the  piezometric  measurement 
wells  in  relation  to  the  simulation  grid. 


Table  2 

Piezometric  Head  Manual  Measurement  Dates 

June  19,  1990 

November  7,  1990 

April  4,  1991 

August  19,  1991 

July  23,  1990 

December  5,  1990 

May  10,  1991 

September  1 1,  1991 

August  13,  1990 

January  8,  1991 

May  20,  1991 

September  1 7,  1990 

February  8,  1991 

June  13,  1991 

October  15,  1990 

March  8,  1991 

July  9,  1991 

Porosity 

The  average  porosity  from  84  core  samples  collected  from  4  coreholes  at 
the  MADE  site  was  0.32  (Boggs  et  al.  1990,  Table  3  below).  These  samples 
were  taken  from  areas  relatively  close  to  the  injection  point.  Three  of  the 
core  sample  locations  were  within  approximately  25  m  of  the  injection  area 
while  the  other  core  sample  location  was  less  than  80  m  from  the  injection 
area. 


Table  3 

Statistical  Summary  of  Porosity  Measurements 

Corehole 

Mean  Porosity 

Standard  Deviation 

CD1 8 

0.30 

0.07 

CD21 

0.32 

0.09 

CM  66 

0.32 

0.09 

CM151 

0.32 

0.07 
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Figure  31 .  Piezometric  well  locations  relative  to  grid 
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Recharge 

The  recharge  rate  for  this  study  was  set  to  a  constant  value  of 
7.0E-04  m/day  using  a  flow  calibration  procedure.  Available 
meteorological  data  from  Memphis,  TN,  Birmingham,  AL,  and  Meridian, 
MS  used  with  the  Environmental  Protection  Agency’s  Pesticide  and  Root 
Zone  Model  (PRZM-2,  Mullins  et  al.  1993)  calculated  an  average  daily 
recharge  rate  of  9.6E-03  m/day  for  those  areas  surrounding  the  site  area. 


Boundary  conditions 

The  lowermost  layer  of  the  grid  was  set  to  a  no-flow  boundary  to  assume 
the  site  conditions  of  a  relatively  impervious  layer.  Boundary  conditions 
used  for  the  328-day  simulation  were  specified-head  boundaries.  The  head 
contours  were  overlain  with  the  MADE  site  grid  and  the  head  for  selected 
boundary  cells  was  obtained  for  the  February  8,  1991  data  set.  The 
February  1991  piezometric  head  data  set  was  chosen  because  it  matched 
closely  with  the  average  head  distribution  of  all  the  data  sets.  The  specified 
head  was  then  placed  into  a  file  representing  the  heads  for  a  selection  of 
outer  boundary  cells  located  on  the  four  vertical  sides  of  the  grid.  A  series 
of  simulations  were  conducted  and  the  specified  boundary  heads  were 
adjusted  until  a  satisfactory  match  between  simulated  and  observed  head  at 
each  measurement  well  location  within  the  simulated  domain  was  obtained. 
In  addition,  a  series  of  simulations  were  run  using  time-variant  boundary 
conditions.  The  results  are  given  in  Appendix  B. 


Stress  periods 

Two  stress  periods  were  defined  in  the  GMS/MODFLOW  model.  This 
corresponds  to  one  stress  period  for  the  injection  period  and  one  stress 
period  for  the  recovery  period.  The  first  stress  period  covered  2  days  of 
injection,  while  the  second  was  326  days  long,  up  to  the  end  of 
measurements. 


MODFLOW  Simulation  Results 


A  series  of  modeling  scenarios  were  run  using  the  defined  conceptual 
model  with  proper  initial  and  boundary  conditions.  The  strongly  implicit 
numerical  solution  option  (SIP)  of  MODFLOW  was  selected  for  all 
applications.  The  simulation  was  performed  using  the  transient  option  of 
MODFLOW.  Two  transient  simulations  were  done.  In  the  first  one, 
time-variable  boundary  conditions  were  applied  and  in  the  second  one, 
time-invariant  average  head  boundary  conditions  were  used.  Results  of  the 
first  approach  were  not  used  in  the  MT3D  simulations  (Chapter  6);  however, 
they  are  given  in  Appendix  B.  The  maximum  number  of  iterations  was  set 
to  600  with  an  acceleration  parameter  of  one.  The  head  change  criterion  for 
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convergence  was  set  to  1.0E-05  m.  The  leakance  required  by  the  model  was 
calculated  using  vertical  hydraulic  conductivities  of  l/10th  the  value  of  the 
horizontal  hydraulic  conductivities.  Storage  coefficients  used  in  the 
simulation  were  assumed  to  be  1.0E-01  for  the  unconfined  model  layer  and 
1.0E-04  for  each  of  the  remaining  layers. 

The  injection  of  chemicals  in  the  field  occurred  through  five  wells 
spaced  1  m  apart  and  screened  over  a  0.6-m  interval  between  elevations  of 

57.5  and  58.1  m.  For  MODFLOW  simulation,  the  five  wells  were  combined 
into  one  well,  since  the  cell  size  of  the  grid  system  covers  the  area  of  the 
five  wells  used  at  the  site.  The  total  fluid  injection  was  9.7  m3  over  a  period 
of  48.5  hr.  This  was  simulated  as  4.85  m3/day  over  a  period  of  2  days. 
Injection  was  induced  into  layer  5  of  the  model,  with  a  top  elevation  of 

58.5  m  and  a  bottom  elevation  of  57.5  m. 

The  MODFLOW  available  in  GMS  Version  1.2  could  not  be  used  to 
simulate  time-variant  specified-head  boundary  conditions.  The 
specified-head  boundary  condition  capability  in  MODFLOW  was  achieved 
by  adding  a  suitable  algorithm  (Leake  and  Prudic  1988)  to  the  MODFLOW 
code.  This  addition  to  the  model  allowed  simulation  of  the  region  of 
interest  with  very  good  results.  A  preprocessor  had  to  be  created  to 
generate  the  file  structure  required  by  this  algorithm  and  is  transparent  to 
GMS. 

The  output  from  a  MODFLOW  simulation  illustrating  the  simulated  head 
for  layers  3  and  7  at  a  time  of  328  days  is  shown  in  Figure  32.  The  figure 
shows  very  steep  gradients  in  the  vicinity  of  the  injection  area  and  the 
leveling  off  of  the  gradients  in  the  mid-field  and  far-field  areas,  as  was 
indicated  from  the  observed  data. 


MODFLOW  Model  Calibration 


As  described  earlier,  hydraulic  conductivity  plays  a  major  role  in  the 
numerical  simulation  of  this  heterogeneous  site.  Discrete  measurements  of 
hydraulic  conductivity  were  interpolated  to  generate  values  at  each  of  the 
grid  blocks.  As  part  of  a  calibration  process,  the  interpolated  hydraulic 
conductivity  and  boundary  conditions  were  changed  manually  within 
reasonable  limits.  Manual  changes  of  these  parameters  continued  until 
differences  between  observed  and  measured  hydraulic  heads  were  not 
significant. 

Table  4  compares  simulated  and  observed  heads  for  layer  2  of  the  model 
grid.  Likewise,  Tables  5  and  6  compare  heads  for  model  layers  3  and  7, 
respectively.  The  maximum  difference  between  simulated  and  observed 
heads  was  0.25  m,  or  0.4  percent,  as  shown  in  Tables  4-6.  The  0.4-percent 
maximum  difference  occurred  at  a  well  located  in  layer  7.  These  layers  were 
selected  for  comparison  because  the  piezometric  measurement  locations 
predominantly  fell  within  these  three  layers  (Figures  33a  and  33b).  It 
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Figure  32.  MODFLOW  simulated  total  head  (m)  for  layers  3  and  7  at  time  =  328  days 


33b).  It  should  be  noted  that  in  this  study,  the  relatively  good  calibrated 
flow  model  was  not  the  best  one  for  the  transport  model.  The  final 
hydraulic  conductivity  profile  used  was  chosen  because  it  demonstrated  the 
best  overall  agreement  between  observed  and  simulated  heads  and  observed 
and  simulated  contaminant  concentrations.  Figure  34  shows  the  February  8, 
1991  observed  piezometric  heads  in  layers  3  and  7  of  the  numerical  grid 
(shallow  and  deep  part  of  aquifer,  respectively).  Although  the  simulated 
results  are  not  shown  in  the  figure,  there  is  a  reasonable  agreement  between 
the  simulated  and  observed  data.  The  cumulative  MODFLOW  volume 
balance  discrepancy  for  the  simulation  was  almost  0  percent. 


Table  4 

Observed  (Obs)  and  Simulated  (Sim)  Heads  in  Model  Layer  2 

Weil  No. 

Obs  (m) 

Sim  (m) 

Difference  (m) 

2 

62.65 

62.78 

+0.13 

3 

62.53 

62.53 

±0.0 

8 

62.62 

62.60 

-0.01 

10 

62.47 

62.54 

+0.07 

17 

62.25 

62.20 

-0.05 

Table  5 

Observed  (Obs)  and  Simulated  (Sim)  Heads  in  Model  Layer  3 

Well  No. 

Obs  (m) 

Sim  (m) 

Difference  (m) 

9 

63.11 

63.03 

-0.08 

11 

62.58 

62.59 

±0.01 

7 

62.52 

62.54 

+0.02 

6 

62.50 

62.53 

+0.03 

13 

62.46 

62.49 

+0.03 

12 

62.49 

62.49 

±0.0 

18 

62.41 

62.44 

+0.03 

15 

62.42 

62.39 

-0.03 

16 

62.37 

62.35 

-0.02 
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Table  6 

Observed  (Obs)  and  Simulated  (Sim)  Heads  in  Model  Layer  7 


Well  No. 

Obs  (m) 

Sim  (m) 

Difference  (m) 

22 

62.80 

62.77 

-0.03 

34 

62.61 

62.64 

+0.03 

26 

62.55 

62.54 

-0.01 

27 

62.50 

62.54 

+0.04 

00 

CO 

62.37 

62.53 

+0.16 

36 

62.46 

62.49 

+0.03 

!  35 

62.49 

62.48 

-0.01 

41 

62.46 

62.44 

-0.02 

;  39 

62.37 

62.35 

-0.02 

!  40 

61.87 

62.12 

+0.25 

MODFLOW  Model  Sensitivity  Analysis 

The  accuracy  of  the  MODFLOW  results  (as  for  any  numerical  model) 
depends  on  the  accuracy  of  measured  or  estimated  input  parameters.  For 
some  applications,  it  is  useful  to  develop  a  criterion  to  estimate  the 
sensitivity  of  the  model  results  to  the  changes  of  input  parameters.  This  will 
provide  a  tool  for  the  user  to  quantify  the  range  of  uncertainties  associated 
with  the  input  parameters.  In  what  follows,  some  of  the  important  input 
parameters  have  been  changed  arbitrarily  and  effects  on  the  simulations  are 
assessed. 


Hydraulic  conductivity 

As  described  earlier,  hydraulic  conductivity  has  a  major  effect  on  the 
outcome  of  simulations.  All  values  of  hydraulic  conductivity  were  increased 
by  10  percent  and  compared  to  simulated  results.  Figure  35  shows  the 
difference  in  simulated  head  when  the  hydraulic  conductivity  values  were 
increased  by  10  percent.  The  hydraulic  head  was  affected  more  in  the  area 
of  converging  flow  and  steep  gradients. 


Leakance 

Leakance  in  MODFLOW  is  a  function  of  the  vertical  hydraulic 
conductivity  and  the  confining  layer  thickness  and  is  needed  for  every  layer 
except  the  bottom  layer  of  the  grid  system.  A  decrease  of  the  leakance  by 
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Figure  33.  Plan  (X-Y)  and  vertical  view  (Z-Y)  of  piezometric  head  measurement  locations 
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Figure  34.  Observed  piezometric  heads  in  layers  3  and  7  of  the  numerical  grid 
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Figure  35.  Differences  in  total  head  in  layers  3  and  7  (time  =  328  days)  caused  by  a  10-percent  increase  in  hydraulic  conductivity 


10  percent  caused  a  change  of  5  percent  in  the  hydraulic  head.  The  changes 
were  noticeable  in  the  area  with  steep  gradients.  This  change  of  leakance 
had  little  effect  on  the  mid-field  and  far-field  regions  of  the  model  grid 
(Figure  36  for  layers  3  and  7)  because  of  horizontal  hydraulic  conductivity 
distribution  in  those  areas. 


Recharge 

Increasing  the  recharge  rate  by  an  arbitrary  factor  of  6.5  changed  the 
head  by  a  maximum  absolute  value  of  less  than  1  percent.  Figures  37,  38, 
and  39  illustrate  the  change  in  head  for  layers  1,  3  and  7,  respectively, 
caused  by  this  increase  in  recharge. 


Summary  of  MODFLOW  Simulations 

The  study  showed  that  one  of  the  most  important  controlling  factors  for 
modeling  a  heterogeneous  system  such  as  CAFB  is  the  numerical 
characterization  of  the  hydraulic  conductivity  field.  A  large  portion  of  the 
time  devoted  to  this  study  was  spent  in  attempting  to  determine  as 
accurately  as  possible  the  hydraulic  conductivity  for  each  cell  of  the  model 
grid.  This  task  was  enormous  in  that  several  commercial  interpolation 
applications  were  used  in  trying  to  determine  the  hydraulic  conductivity  to 
be  used  in  the  model  grid.  In  addition  to  the  various  interpolation  packages 
used,  several  interpolation  methods  were  attempted.  Another  controlling 
factor  in  the  simulation  is  the  grid  size  in  relation  to  the  variance  of  the 
hydraulic  conductivity.  As  mentioned  previously,  the  grid  size  was  5m  x 
5m  x  lm  and  the  variance  of  the  natural  logarithm  of  hydraulic  conductivity 
was  4.5.  The  extreme  heterogeneity  of  the  site  requires  a  finer  resolution 
grid  to  achieve  better  results.  However,  a  finer  resolution  grid  in  turn 
requires  additional  computational  effort  and  increased  storage  capacity  for 
the  simulations.  A  nonuniform  grid  perhaps  could  provide  better  results 
with  less  computer  storage  requirements.  However,  if  a  detailed  correct 
distribution  of  hydraulic  conductivity  could  not  be  defined,  the  nonuniform 
grid  would  not  be  better  than  the  uniform  grid.  The  accuracy  desired, 
computational  effort,  and  storage  requirements  are  factors  to  consider  in  any 
particular  application  of  GMS/MODFLOW. 

The  flow  simulation  appeared  to  match  the  observed  data  relatively  well. 
The  simulated  piezometric  heads  in  most  cases  were  within  0.2  m  of  the 
observed  data.  The  largest  difference  obtained  between  observed  and 
simulated  values  was  0.25  m  (less  than  1  percent).  The  simulation  assumed 
steady-state  conditions  taken  as  the  average  of  the  piezometric  head 
monitoring  data  from  June  1990  through  September  1991.  The  average 
piezometric  head  for  this  time  period  was  very  close  to  the  February  1991 
monitoring  data  and  thus  these  data  were  used  to  drive  the  flow  model.  The 
water  table  at  the  site  had  about  20  to  30  percent  variation  between  low-  and 
high-water-table  seasons.  However,  the  simulations  conducted  assumed  an 
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Figure  36.  Change  in  total  head  (m)  in  layers  3  and  7  due  to  leakance  change 
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Figure  37.  Change  in  total  head  (m)  in  layer  1  due  to  increase  of  recharge  rate 
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Figure  38.  Changes  in  total  head  (m)  in  layer  3  due  to  increase  of  recharge  rate 
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Figure  39.  Change  in  total  head  (m)  in  layer  7  due  to  increase  of  recharge  rate 


an  average  constant  head  for  the  duration  of  the  simulation.  It  is  very  likely 
that  the  simulation  could  have  been  improved  overall  by  simulating  a  low 
and  high  water  table  for  the  simulation  period  of  328  days  as  described  in 
the  FEMWATER  application  chapters.  The  simulation  without  the 
high-water-table  event  showed  that  the  travel  extent  of  simulated 
contaminant  plumes  is  limited  compared  to  observed  plumes. 

The  simulation  did  not  appear  to  be  very  sensitive  to  the  recharge  rate 
from  the  scenarios  evaluated.  Increasing  the  rate  of  recharge  by  6.5  times 
caused  a  maximum  absolute  change  in  head  of  less  than  1  percent. 

The  model  was  sensitive  to  the  changes  in  hydraulic  conductivity  and 
indicates  the  importance  of  a  good  hydraulic  conductivity  representation  for 
the  model.  Increasing  the  hydraulic  conductivity  for  all  cells  by  10  percent 
caused  a  notable  difference  in  the  simulated  head  in  the  area  near  the 
injection  point.  The  effects  were  almost  totally  limited  to  the  near-field 
region  of  the  simulation  domain.  There  was  little  effect  of  the  hydraulic 
conductivity  change  in  the  middle-  and  far-field  regions. 

The  model  was  also  sensitive  to  changes  of  leakance.  Decreasing  the 
leakance  by  10  percent  caused  a  maximum  absolute  change  of  5  percent  in 
the  pressure  head.  The  effect  of  this  change  was  most  notable  in  the  lower 
layers  and  also  in  the  far-field  region  of  the  simulated  domain.  The 
far-field  region  coincides  with  steep  vertical  gradients. 

Calibration  was  done  manually  by  changing  the  hydraulic  conductivity 
and  boundary  condition  pressure  heads.  This  could  be  better  accomplished 
by  the  use  of  a  parameter  estimation  program.  The  program  should  be 
capable  of  running  the  simulation,  comparing  the  simulated  values  to 
observed  values,  and  then  revising  the  specified-head  boundary  condition 
file  and  repeating  this  process  until  a  satisfactory  match  can  be  obtained.  In 
regard  to  the  MODFLOW  model,  MODFLOWP  (Hill  1990)  provides  these 
capabilities. 

Recommendations  for  future  simulations  include  better  numerical 
characterization  of  the  hydraulic  conductivity  for  the  site  and  a  model  grid 
with  finer  resolution  and/or  an  adaptive  grid.  Models  capable  of  simulating 
both  the  saturated  and  unsaturated  zones  with  transient  boundary  conditions 
perhaps  provide  better  simulation  of  the  site  flow  conditions. 
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6  MT3D  Modeling  Approach 


This  chapter  describes  the  application  of  a  modular  three-dimensional 
transport  model,  MT3D/GMS  (Zheng  1990),  to  the  CAFB  site.  The  selected 
chemicals  used  in  the  simulations  are  tritium,  14C,  benzene,  p-xylene, 
naphthalene,  and  o-dichlorobenzene.  The  chemical  characteristics  of  the 
groundwater  at  the  site  are  described  in  Boggs  et  al.  (1990).  The  flow 
information  required  by  MT3D  was  calculated  previously  using 
MODFLOW  as  described  in  Chapter  5. 


Summary  of  MT3D  Model 

MT3D  is  a  3-D  solute  transport  model  which  calculates  advection, 
dispersion,  and  chemical  reactions  (1st  order  decay)  of  contaminants  in 
groundwater  systems.  The  numerical  technique  used  in  MT3D  is  a 
block-centered  finite  difference.  The  MT3D  computer  program  adopts  a 
modular  structure  to  permit  simulation  of  transport  components 
independently  or  jointly.  MT3D  is  designed  for  linkage  to  the  MODFLOW 
model  (Chapter  5).  However,  other  flow  models  also  may  be  adopted  for 
use  by  MT3D.  MT3D  has  three  basic  solution  options:  the  method  of 
characteristics  (MOC),  modified  method  of  characteristics  (MMOC),  and  a 
hybrid  of  these  two  methods  (HMOC).  This  approach  combines  the  strength 
of  MOC  for  reducing  numerical  dispersion  and  the  computational  efficiency 
of  MMOC,  making  MT3D  flexible  for  a  wide  range  of  field  problems.  In 
addition,  MT3D  can  also  be  used  in  a  standard  finite-difference  mode.  The 
general  types  of  boundary  conditions  can  be  assigned  in  MT3D.  The 
boundary  conditions  can  be  either  a  known  concentration  around  the 
boundary  (Dirichlet),  a  concentration  gradient  across  the  boundary 
(Neumann),  or  a  combination  of  these  two. 
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MT3D  Model  Input 


Flow  data 

MT3D/GMS  requires  inter-cell  flow/flux  data  calculated  by  MODFLOW. 
The  required  flow  data  are  stored  in  a  linkage  file  that  is  assigned  by  the 
user  during  the  MODFLOW  application.  The  flow  parameters  used  to  drive 
the  transport  model  are  described  in  Chapter  5. 


Chemical  injection  data 

A  total  of  9.7  m3  of  the  selected  chemical  solution  was  injected  into  the 
aquifer  over  a  2-day  period.  The  injection  concentrations  and  total  mass 
injected  are  summarized  in  Table  7.  The  tracer  solution  was  prepared  from 
ambient  groundwater  taken  75  m  up-gradient  from  the  injection  point.  The 
tritium  and  pure-phase  organic  compounds  were  mixed  and  stored  in  two 
2.7-m3  thermally  insulated  tanks  equipped  with  a  propeller  mixer  and 
floating  lids.  Samples  of  the  tracer  solution  were  collected  at  1-  and  2-hr 
intervals  from  tubes  leading  to  the  injection  wells.  The  mass  injected  for 
each  tracer  was  estimated  by  volume-weighting  of  these  tracer 
concentration  samples. 


Table  7 

Injected  Chemicals 

Tracer 

Mean  Concentration 

Mass  Injected 

tritium 

55,610  pCi/ml 

0.5387  Ci 

14C  (p-xylene) 

2,770  pCi/ml 

0.0268  Ci 

benzene 

68.1  mg/L 

659.7  g 

p-xylene 

51.5  mg/L 

402.0  g 

o-dichlorobenzene 

32.8  mg/L 

317.7  g 

naphthalene 

7.23  mg/L 

70.0  g 

The  downstream  transport  of  these  contaminants  was  tracked  by 
sampling  from  multiple  ports  in  wells  for  a  given  time  period  using  an  array 
of  328  multi-level  samplers.  Details  of  the  sampling  procedure  are  given  in 
a  report  by  EPRI  (1993b).  The  samplings  were  done  27,  132,  224,  and  328 
days  after  injection.  A  special  sampling  strategy  was  developed  for  the 
organic  tracer  since  volatile  losses  from  the  former  sampling  method  could 
reduce  organic  tracer  concentrations.  This  alternative  sampling  method 
used  56  positive-displacement  samplers  oriented  along  two  parallel  rows  or 
“fencelines”  normal  to  the  flow  direction  in  the  area  near  the  injection  site. 
The  rows  were  placed  10  m  apart  and  were  located  at  distances  of 
approximately  6  and  16  m  from  the  tracer  injection  point.  The  sampling 
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devices  were  installed  in  15-cm-diam  auger  holes  with  no  more  than  two 
sampling  points  per  hole.  The  sampling  points  were  isolated  by  bentonite 
seals  located  approximately  0.6  m  above  and  below  the  samplers  in  the 
annular  space.  The  sampling  frequency  for  the  fencelines  varied  from 
2-week  intervals  to  3-month  intervals. 

Boggs  et  al.  (1993)  estimated  the  biodegradation  rates  from  mass  balance 
analysis  of  MADE-2  data  (Table  8).  These  data  were  used  as  the  first-order 
decay  coefficient  in  the  MT3D  model. 


Table  8 

Biodegradation  Rates  Used  in  MT3D  Simulations  | 

T  racer 

Estimated  rate  (d'1) 

Lower  95% 
confidence  interval 

Upper  95% 
confidence  interval 

benzene 

0.008 

0.004 

0.011 

p-xylene 

0.011 

0.006 

0.016 

o-dichlorobenzene 

0.006 

0.001 

0.010 

naphthalene 

0.007 

0.004 

0.010 

The  sorption  data  used  in  this  study  were  obtained  from  laboratory  data 
obtained  from  prior  studies  on  the  MADE  site  (Boggs  et  al.  1993).  Table  9 
indicates  the  sorption  data  used  in  this  study  for  benzene,  naphthalene, 
p-xylene,  and  o-dichlorobenzene. 


Table  9 

Sorption  Data  Used  in  MT3D  Simulations 

T  racer 

Kd  (mL/g) 

benzene 

0.059 

p-xylene 

0.048 

o-dichlorobenzene 

0.065 

naphthalene 

0.085 

Model  grid 

The  uniform  grid  developed  with  GMS  for  the  MADE  site  contaminant 
transport  study  consisted  of  25,025  block  elements  in  11  layers,  65  rows, 
and  35  columns,  as  shown  in  Figure  29.  The  elevation  of  the  top  of  the  first 
layer  was  set  to  62.5  m  with  the  thickness  of  each  underlying  layer  being  1 
m  and  a  bottom  elevation  of  the  lowermost  layer  of  51.5  m. 
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Boundary  conditions 


The  lowermost  layer  of  the  grid  coincided  with  the  no-flow  boundary  of 
the  flow  grid  and  reflects  the  site  conditions  of  an  impervious  layer 
approximately  at  this  location.  The  boundary  conditions  used  for  the 
328-day  simulation  were  no-mass  flux  boundaries. 


Stress  periods 

The  simulation  used  two  stress  periods.  This  corresponds  to  one  stress 
period  for  the  injection  and  one  stress  period  for  the  recovery.  The  total 
simulation  time  for  the  model  was  328  days. 


Flow  linkage 

The  MODFLOW  results  were  required  as  input  to  the  MT3D  transport 
model.  The  actual  injection  in  the  field  occurred  through  five  wells  spaced 
1  m  apart  and  screened  over  a  0.6-m  interval  between  elevations  of  57.5  and 
58.1  m.  The  actual  field  injection  rate  for  the  five  wells  was  3.3  L/min  over 
a  period  of  48.5  hr  with  a  total  volume  injected  of  9.7  m3.  For  the 
simulations,  the  flow  rates  of  five  wells  in  the  field  were  combined  into  a 
single  well  with  a  flow  rate  of  4.8  m3/day.  Injection  occurred  in  layer  5  of 
the  model,  which  had  a  top  elevation  of  58.5  m  and  a  bottom  elevation  of 
57.5  m. 


Transport  simulations 

Transport  simulations  incorporated 
(a)  advection  and  dispersion,  and  (b) 
advection-only  options  of  MT3D  into  the 
calculations.  Advection  was  attempted  with 
several  solution  scheme  options  including 
MOC,  MMOC,  HMOC,  and  upstream  finite 
difference.  It  appeared  that  the  hybrid  method 
of  characteristics  advective  solution  scheme 
and  no  dispersion  may  provide  better  results 
for  this  site.  The  particle  tracking  algorithm 
used  for  the  study  was  Runge-Kutta  at  or  near 
sources  and  Euler  elsewhere  with  a 
concentration  weighting  factor  of  0.5.  The 
number  of  particles  and  the  maximum  number 
of  cells  a  particle  was  allowed  to  move  during 
the  simulations  are  shown  in  Figure  40. 
Porosity  for  the  simulation  grid  was  set  to  0.10 
in  the  meander  channel  and  0.31  elsewhere. 
The  simulation  of  each  of  the  five  modeled 
contaminants  is  described  below. 
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Figure  40.  View  of  GMS/MT3D  particle 
screen  showing  parameters 
used  in  simulations 
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Simulation  of  tritium 


MT3D  was  used  to  simulate  328  days  of  tritium  transport  at  the  CAFB 
site.  Tritium  is  treated  as  a  conservative  tracer,  i.e.,  negligible  adsorption 
and  decay. 

Output  from  MT3D  simulation  expressed  in  plan  view  and  in  oblique 
view  (vertical  view)  for  a  50-pCi/mL  iso-surface  is  shown  in  Figures  41  and 
42,  respectively.  Figure  43  shows  an  oblique  view  of  a  5-pCi/ml  tritium 
iso-surface  at  days  27,  132,  224,  and  328.  The  observed  concentrations  of 
tritium  are  shown  in  Figures  44  and  45  (  Boggs  et  al.  1993).  Note  that 
Figures  44  and  45  show  the  relative  concentration,  and  the  scale  of  these 
graphs  is  different  from  the  one  used  in  Figures  41,  42,  and  43. 

When  the  simulated  and  observed  plumes  are  compared  graphically, 
several  observations  as  to  the  relative  agreement  between  the  plumes  can  be 
made.  At  day  27,  the  simulated  and  observed  plumes  agree  well  with  each 
other,  and  velocity  and  low  hydraulic  conductivity  in  the  region  are 
mimicked  well.  The  simulated  and  observed  plumes  of  50-pCi/mL 
iso-surface  at  day  132  agree  well  with  each  other.  However,  the  observed 
and  simulated  plumes  at  a  lower  concentration  of  5  pCi/mL  showed  a  larger 
discrepancy  in  terms  of  the  areal  extent  of  the  plume  (the  simulated  lagging 
approximately  56  m  behind  the  observed).  At  day  224,  the  observed 
5-pCi/mL  plume  had  reached  a  downgradient  extent  in  excess  of  270  m 
while  the  simulated  plume  had  only  traveled  approximately  80  m  from  the 
initial  point  of  injection.  The  50-pCi/mL  observed  and  simulated  plumes 
have  traveled  nearly  the  same  extent.  The  observed  5-pCi/mL  and  the 
50-pCi/mL  iso-surfaces  at  this  time  period  are  noticeably  different  in  travel 
extent.  The  observed  plume  differences  suggest  the  possibility  of 
preferential  flow  pathways  that  the  monitoring  equipment  could  not  discern 
due  to  monitoring  location.  At  day  328,  the  simulated  plume  of  5-pCi/mL 
iso-surfaces  lagged  the  observed  plume  by  approximately  78  m.  The 
50-pCi/mL  simulated  and  observed  iso-surfaces  show  the  simulated  plume 
having  traveled  approximately  177  m  further  than  the  observed  plume. 

Simulation  of  carbon  14  (14C) 

The  simulation  of  14C  transport  through  the  MADE  system  was 
performed  for  a  period  of  328  days.  The  decay  rate  used  for  carbon-14  was 
based  on  a  half-life  of  5,700  years  and  no  adsorption  was  simulated  for  this 
tracer.  Output  from  the  MT3D  simulation  expressed  in  plan  view  for  a 
3-pCi/mL  iso-surface  is  shown  in  Figure  46  for  days  27,  132,  242,  and  328. 
The  simulated  3-pCi/mL  iso-surface  is  shown  in  oblique  view  in  Figure  47 
for  days  27,  132,  234,  and  328.  Figure  48  shows  the  simulated  plume  with  a 
concentration  of  0.3  pCi/mL.  Figure  49  (Boggs  et  al.  1993)  shows  the 
vertical  and  longitudinal  extent  of  the  relative  concentration  14C  plume  at 
27,  132,  224,  and  328  days.  The  14C  plume  and  tritium  plumes  show  very 
similar  behavior,  which  was  expected  and  is  evident  in  the  observed  plume 
figures.  The  14C  plume  does  not  travel  the  downgradient  extent  of  the 
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Figure  41. 


Plan  view  of  simulated  50-pCi/ml  tritium  iso-surface  at  a)  day  27,  b)  day  132,  c)  day  234,  and 
d)  day  328 
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Figure  43.  Oblique  view  of  simulated  5-pCi/ml  tritium  iso-surface  at  a)  day  27,  b)  day  132,  c)  day  234,  and  d)  day  328 
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Figure  45.  Observed  vertical  and  longitudinal  extent  of  tritium  plume  at  27,  132,  224,  and  328  days 
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Figure  46.  Plan  view  of  simulated  3-pCi/ml  carbon-14  iso-surface  at  a)  day  27, 
b)  day  132,  c)  day  234,  and  d)  day  328 
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observed  plume  nor  does  it  show  the  spreading  of  the  advancing  front.  As 
shown  in  the  figures  of  the  simulated  and  observed  14C  plumes  at  27,  132, 
224,  and  328  days,  the  simulation  was  able  to  match  reasonably  well  until 
day  224,  when  the  simulation  lagged  the  observed  plume.  At  day  328  the 
simulation  of  the  0.3-pCi/ml  iso-surface  still  lagged  the  observed  plume, 
but  showed  a  tendency  to  move  further  downgradient,  probably  due  to 
increased  velocity  in  this  area  of  the  aquifer.  Also  at  day  328,  the  simulated 
and  observed  3-pCi/ml  plumes  showed  that  the  simulated  plume  had 
traveled  further  downgradient  than  the  observed  plume.  The  reason  for  the 
lag  between  the  simulated  plume  and  the  observed  plume  could  be  due  to 
the  hydrologic  conditions  occurring  at  about  this  time.  Piezometric  data 
show  that  at  approximately  day  203  of  the  simulation,  the  piezometric  head 
is  at  its  highest  point  since  the  beginning  of  the  monitoring.  Only 
steady-state  flow  was  conducted  here  and  thus  the  flow  model  lacks  the 
seasonal  fluctuation  in  the  water  table.  The  lower  concentration 
0.3-pCi/mL  observed  and  simulated  plumes  show  a  larger  discrepancy 
(approximately  56  m)  between  them.  At  day  224,  the  simulated  and 
observed  0.3-pCi/mL  iso-surfaces  have  significant  differences.  The 
observed  0.3-pCi/mL  plume  has  reached  a  downgradient  extent  in  excess  of 
270  m  from  the  injection  point,  while  the  simulated  plume  has  traveled  only 
approximately  80  m  from  the  initial  point  of  injection.  The  3-pCi/mL 
observed  and  simulated  plumes  have  traveled  nearly  the  same  extent.  The 
difference  between  the  observed  0.3-pCi/mL  and  the  3-pCi/mL  iso-surfaces 
at  this  time  period  is  noticeably  different  in  travel  extent.  The  observed 
plume  differences  suggest  the  possibility  of  preferential  flow  pathways  that 
the  monitoring  equipment  could  not  discern  due  to  monitoring  location.  At 
day  328,  the  observed  and  simulated  0.3-pCi/mL  iso-surfaces  show  the 
simulated  plume  lagging  the  observed  plume  by  approximately  78  m.  The 
3-pCi/mL  simulated  and  observed  iso-surfaces  show  the  simulated  plume 
having  traveled  approximately  177  m  further  than  the  observed  plume. 


Simulation  of  benzene 

Benzene  transport  through  the  MADE  system  was  simulated  for  a  period 
of  328  days.  Input  parameters  included  adsorption  and  first-order  biodecay 
coefficients  (Tables  7,  8,  and  9).  Output  from  the  MT3D  simulation 
expressed  in  plan  view  for  a  50-mg/m3  iso-surface  is  shown  in  Figure  50  for 
days  27,  132,  242,  and  328.  The  simulated  50-mg/m3  iso-surface  is  shown  in 
oblique  view  in  Figure  51  for  days  27,  132,  234,  and  328.  Figure  52  shows 
an  oblique  view  of  a  5-mg/m3  benzene  iso-surface  at  days  27,  132,  224,  and 
328.  Figure  53  shows  the  observed  vertical  and  longitudinal  extent  of  the 
relative  concentration  benzene  plume  at  27,  132,  224,  and  328  days  (Boggs 
et  al.  1993).  The  simulated  benzene  plume  advanced  downgradient  more 
than  the  observed  plume,  but  overall  was  a  fairly  good  match.  As  can  be 
seen  from  the  figures  of  the  simulated  and  observed  benzene  plumes  at  27, 
132,  224,  and  328  days,  the  simulation  was  able  to  match  reasonably  well 
until  day  224,  after  which  the  simulated  plume  propagated  further  than  the 
observed  plume.  At  day  328,  the  simulated  5-mg/m3  and  50-mg/m3 
iso-surfaces  had  still  traveled  further  than  the  observed  corresponding 
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Figure  50.  Plan  view  of  simulated  50-mg/m  benzene  iso-surface  at  a)  day  27 
b)  day  132,  c)  day  234,  and  d)  day  328 
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Figure  53.  Observed  vertical  and  longitudinal  extent  of  benzene  plume  at  27, 
132,  224,  and  328  days 
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plumes,  but  had  shown  a  tendency  to  remain  in  place  since  it  had  only 
moved  approximately  10  m  further  downgradient  over  a  104-day  period.  At 
day  27,  the  simulated  and  observed  plumes  agree  well  with  each  other.  The 
simulated  and  observed  50-mg/m3  plumes  both  traveled  downgradient 
approximately  25  m  while  the  5-mg/m3  iso-surfaces  had  both  traveled 
approximately  25  m  also.  At  day  132,  the  simulated  and  observed  plumes 
of  the  50-mg/m3  and  5-mg/m3  iso-surfaces  agree  well  with  each  other.  At 
day  224,  the  simulated  and  observed  50-mg/m3  and  5-mg/m3  iso-surfaces 
depart  from  the  behavior  of  the  tritium  and  14C  contaminants  in  that  the 
simulated  plumes  travel  further  downgradient  than  the  observed  plumes. 

The  observed  50-mg/m3  plume  has  reached  a  downgradient  extent  of 
approximately  15  m  while  the  simulated  plume  traveled  approximately  40  m 
from  the  initial  point  of  injection.  At  day  328,  the  observed  50-mg/m3 
iso-surface  shows  the  observed  plume  has  still  remained  at  approximately 
15  m  from  the  injection  point,  while  the  simulated  plume  has  extended  to 
approximately  75  m  from  the  source.  The  5-mg/m3  simulated  and  observed 
iso-surfaces  show  the  simulated  plume  having  traveled  approximately  70  m 
further  than  the  observed  plume. 


Simulation  of  p-xylene 

P-xylene  transport  in  the  CAFB  site  was  simulated  for  a  period  of  328 
days.  The  input  parameters  include  adsorption  and  biodecay  coefficients 
(Tables  7,  8,  and  9).  Output  from  the  MT3D  simulation  expressed  in  plan 
view  for  a  50-mg/m3  iso-surface  is  shown  in  Figure  54  for  days  27,  132, 

242,  and  328.  The  simulated  50-mg/m3  iso-surface  is  shown  in  oblique  view 
in  Figure  55  for  days  27,  132,  234,  and  328.  Figure  56  shows  an  oblique 
view  of  a  5-mg/m3  p-xylene  iso-surface  at  days  27,  132,  224,  and  328. 
Figure  57  (Boggs  et  al.  1993)  shows  the  vertical  and  longitudinal  extent  of 
the  relative  concentration  p-xylene  plume  at  27,  132,  224,  and  328  days. 

The  simulated  results  match  reasonably  well  with  observed  data  until  day 
224,  after  which  the  observation  lagged  the  simulated  plume.  At  day  328, 
both  the  simulated  5-mg/m3  and  50-mg/m3  iso-surfaces  had  still  traveled 
further  than  the  observed  corresponding  plumes.  Piezometric  data  show 
that  at  approximately  day  203  of  the  simulation,  the  piezometric  head  is  at 
its  highest  since  the  beginning  of  the  monitoring  and  because  only 
steady-state  flow  was  conducted,  the  flow  model  lacks  the  seasonal 
fluctuation  in  the  water  table.  At  day  27,  the  simulated  and  observed  plumes 
agree  well  with  each  other.  The  simulated  and  observed  50-mg/m3  plumes 
have  a  longitudinal  extent  of  approximately  20  and  25  m  from  the  injection 
point,  respectively.  The  observed  5-mg/m3  plume  has  traveled 
approximately  25  m  from  the  source,  while  the  simulated  plume  has  also 
traveled  approximately  25  m.  At  day  132,  the  simulated  and  observed 
plumes  agree  well  with  each  other  when  comparing  the  50-mg/m3  and 
5-mg/m3  iso-surfaces.  At  day  224,  the  simulated  and  observed  50-mg/m3 
and  5-mg/m3  iso-surfaces  depart  from  the  behavior  of  the  tritium  and  14C 
plumes;  the  simulated  plumes  traveled  further  downgradient  than  the 
observed  plumes.  The  limited  displacement  of  the  observed  plumes  is 
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Figure  54.  Plan  view  of  simulated  50-mg/m3  p-xylene  iso-surface  at  a)  day  27,  b)  day  132,  c)  day  234, 
and  d)  day  328 
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consistent  with  the  effects  of  bio-transformation,  which  prior  studies  at  the 
MADE  site  indicate  are  quite  significant.  The  observed  50-mg/m3  plume 
has  reached  a  downgradient  extent  of  approximately  15  m  while  the 
simulated  plume  traveled  approximately  40  m  from  the  initial  point  of 
injection.  At  day  328,  the  observed  50-mg/m3  iso-surface  showed  the 
observed  plume  has  still  remained  at  approximately  15  m  from  the  injection 
point,  while  the  simulated  plume  extended  to  approximately  75  m  from  the 
source.  The  5-mg/m3  simulated  and  observed  iso-surfaces  show  the 
simulated  plume  traveled  approximately  70  m  further  than  the  observed 
plume. 


Simulation  of  o-dichlorobenzene 

O-dichlorobenzene  transport  through  the  MADE  system  was  simulated 
for  a  period  of  328  days.  Input  parameters  included  adsorption  and 
biodegradation  coefficients  (Tables  7,  8,  and  9).  Output  from  the  MT3D 
simulation  expressed  in  plan  view  for  a  50-mg/m3  iso-surface  is  shown  in 
Figure  58  for  days  27,  132,  242,  and  328.  The  simulated  50-mg/m3 
iso-surface  is  shown  in  oblique  view  in  Figure  59  for  days  27,  132,  234,  and 
328.  Figure  60  shows  an  oblique  view  of  the  5-mg/m3  o-DCB  iso-surface  at 
days  27,  132,  224,  and  328.  Figure  61  shows  the  vertical  and  longitudinal 
extent  of  the  observed  plumes  of  o-dichlorobenzene  at  27,  132,  224,  and 
328  days  (Boggs  et  al.  1993).  Note  that  the  observed  plumes  are  given  as 
nondimensional  relative  concentration  (C/CQ). 

The  simulated  o-dichlorobenzene  plume  showed  behavior  similar  to  that 
of  benzene.  This  behavior  contradicts  the  observed  data  (Figure  61),  which 
shows  that  the  o-DCB  plume  is  similar  to  the  tritium  plume  for  the  first  132 
days.  After  the  first  132  days,  the  observed  plume  deviates  from  the  tritium 
behavior  and  shows  spatial  discontinuity.  As  can  be  seen  from  the  figures 
of  the  simulated  and  observed  o-dichlorobenzene  plumes  at  27,  132,  224, 
and  328  days  (Figures  59  and  61),  the  simulation  was  not  able  to  match  the 
observed  data  for  any  period  beyond  the  first  27  days.  Piezometric  data 
show  that  at  approximately  day  203  of  the  simulation,  the  piezometric  head 
is  at  its  highest  since  the  beginning  of  the  monitoring.  Because  only 
steady-state  flow  was  conducted,  the  flow  input  to  the  model  lacks  the 
seasonal  fluctuation  in  the  water  table  required  to  get  a  better  simulation. 
When  the  simulated  and  observed  plumes  are  compared,  several 
observations  as  to  the  relative  agreement  between  the  plumes  can  be  made. 
At  day  27,  the  simulated  and  observed  plumes  agree  well  with  each  other. 
The  observed  50-mg/m3  plume  has  progressed  downgradient  approximately 
25  m  from  the  injection  point.  The  simulated  50-mg/m3  plume  has  also 
moved  approximately  25  m  from  the  injection  point.  At  day  132,  the 
observed  plumes  progressed  further  downgradient  than  the  simulated 
plumes  when  comparing  the  50-mg/m3  and  5-mg/m3  iso-surfaces.  The 
simulated  5-mg/m3  plume  has  moved  approximately  28  m  from  the  injection 
point,  while  the  observed  plume  has  moved  approximately  65  m.  At  day 
224,  the  simulated  and  observed  50-mg/m3  and  5-mg/m3  iso-surfaces  are 
considerably  different.  The  50-mg/m3  observed  and  simulated  plumes  have 
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Figure  58.  Plan  view  of  simulated  50-mg/m3  o-DCB  iso-surface  at  a)  day  27,  b)  day  132,  c)  day  234,  and 
d)  day  328 
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reached  a  downgradient  extent  of  approximately  60  m  from  the  initial  point 
of  injection. 

The  simulated  5-mg/m3  iso-surface  has  moved  approximately  70  m  from 
the  injection  point,  while  the  observed  plume  has  moved  downgradient  in 
excess  of  260  m  from  the  injection  point.  At  day  328,  the  observed 
50-mg/m3  iso-surface  shows  the  observed  plume  at  approximately  40  m 
from  the  injection  point,  while  the  simulated  plume  has  extended  to 
approximately  75  m  from  the  source.  The  5-mg/m3  simulated  and  observed 
iso-surfaces  show  the  simulated  plume  having  traveled  approximately  85  m 
from  the  injection  point,  while  the  observed  plume  has  traveled 
approximately  130  m  from  the  injection  point. 


Simulation  of  naphthalene 

Naphthalene  transport  through  the  MADE  system  was  simulated  for  a 
period  of  328  days.  Input  parameters  included  adsorption  and  biodecay 
coefficients.  Output  from  the  MT3D  simulation  was  expressed  in  graphs  of 
concentration.  A  plan  view  for  a  7-mg/m3  iso-surface  is  shown  in  Figure  62 
for  days  27,  132,  242,  and  328.  The  simulated  50-mg/m3  iso-surface  is 
shown  in  an  oblique  view  in  Figure  63  for  days  27,  132,  234,  and  328. 
Figure  64  shows  the  oblique  view  of  the  0. 7-mg/m3  naphthalene  iso-surface 
at  days  27,  132,  224,  and  328,  respectively.  The  measured  normalized 
concentration  of  naphthalene  (Boggs  et  al.  1993)  in  the  vertical  and 
longitudinal  extent  at  27,  132,  224,  and  328  days  is  shown  in  Figure  65.  The 
simulated  naphthalene  plume  showed  behavior  similar  to  that  of  the 
observed  plume,  but  again  fails  to  reproduce  the  downgradient  spreading  in 
the  observed  data. 

The  simulated  naphthalene  plumes  matched  the  observed  concentration 
plume  reasonably  well  until  day  224,  after  which  the  observation  lagged  the 
simulation.  At  day  328,  the  simulated  7-mg/m3  and  0. 7-mg/m3  iso-surfaces 
had  still  traveled  further  than  the  observed  corresponding  plumes,  but  had 
shown  a  tendency  to  remain  in  place  since  that  time  by  moving 
approximately  10  m  further  downgradient  over  a  period  of  104  days.  At  day 
27,  the  simulated  and  observed  plumes  agree  well  with  each  other.  The 
observed  7-mg/m3  plume  has  progressed  downgradient  approximately  25  m 
from  the  injection  point. 

The  simulated  7-mg/m3  plume  has  moved  approximately  20  m  from  the 
injection  point.  The  simulated  0. 7-mg/m3  plume  has  moved  approximately 
25  m  from  the  injection  point,  while  the  observed  0. 7-mg/m3  plume  has 
moved  approximately  25  m.  At  day  132,  the  simulated  plumes  have 
progressed  further  downgradient  than  the  observed  plumes  when  comparing 
the  7-mg/m3  and  0. 7-mg/m3  iso-surfaces.  The  simulated  7-mg/m3  plume 
has  moved  approximately  28  m  from  the  injection  point  while  the  observed 
plume  has  moved  approximately  15  m.  At  day  224,  the  simulated  and 
observed  7-mg/m3  and  0. 7-mg/m3  iso-surfaces  are  slightly  different.  The 
simulated  7-mg/m3  plume  has  reached  a  downgradient  extent  of 
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Figure  62.  Plan  view  of  simulated  7-mg/m3  naphthalene  iso-surface  at  a)  day  27,  b)  day  132,  c)  day  234, 
and  d)  day  328 
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Figure  63.  Oblique  view  of  simulated  7-mg/m3  naphthalene  iso-surface  at  a)  day  27,  b)  day  132,  c)  day  234,  and  d)  day  328 
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Figure  65.  Observed  vertical  and  longitudinal  extent  of  naphthalene  plume  at  27, 
132,  224,  and  328  days 
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approximately  50  m  from  the  initial  point  of  injection,  while  the  observed 
7-mg/m3  plume  has  remained  at  approximately  15  m  from  the  source.  The 
simulated  0. 7-mg/m3  iso-surface  has  moved  approximately  60  m  from  the 
injection  point,  while  the  observed  plume  has  moved  downgradient 
approximately  25  m  from  the  injection  point.  At  day  328,  the  observed 
7-mg/m3  iso-surface  shows  the  observed  plume  has  moved  approximately 
20  m  from  the  injection  point,  while  the  simulated  plume  has  extended  to 
approximately  65  m  from  the  source.  The  0. 7-mg/m3  simulated  and 
observed  iso-surfaces  show  the  simulated  plume  having  traveled 
approximately  75  m  from  the  injection  point,  while  the  observed  plume  has 
traveled  approximately  20  m  from  the  injection  point. 

Figure  66  shows  the  GMS  input  screen  for  point  sources  and  shows  the 
injection  concentration  (7,230  mg/m3)  applied  for  the  simulation  of 
naphthalene.  A  linear  adsorption  was  assumed  for  naphthalene.  A  view  of 
the  GMS/MT3D  chemical  reaction  package  is  shown  in  Figure  67. 


Figure  66.  View  of  GMS/MT3D  point  source  entry 
screen 
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MT3D  Model  Calibration 


Generally,  the  calibration  of  the  flow  model  has  been  sufficient  and 
provided  reasonable  input  for  the  transport  model.  In  the  flow  model 
(Chapter  5),  the  model  was  calibrated  for  total  head  distribution.  In  the 
transport  model,  some  of  the  input  parameters  were  changed  in  order  to 
reproduce  observed  plume  behavior. 

Advection  and  dispersion  were  two  main  causes  of  the  plume  spreading 
at  the  site.  The  initial  phase  of  the  simulation  process  used  advection  and 
dispersion.  When  dispersion  was  simulated,  the  contaminant  plumes 
dispersed  too  quickly.  Another  observation  was  that  the  plume  moved 
backward  from  the  natural  gradient  in  such  a  manner  that  it  would  leave  the 
domain  of  the  simulation  grid.  Thus,  modeling  was  conducted  using 
advection  alone  because  it  was  believed  that  the  spatial  variability  of  the 
hydraulic  conductivity  of  the  study  site  would  compensate.  When 
dispersion  was  eliminated,  the  problem  with  the  spreading  of  the  plume 
against  the  natural  gradient  and  outside  the  domain  of  the  simulation  grid 
was  avoided.  Advection  alone  was  able  to  reproduce  reasonably  well  the 
general  shape  of  the  tritium  contaminant  plume  up  until  approximately  day 
224.  This  day  roughly  coincides  with  the  increased  seasonal  fluctuation  of 
the  water  table  (see  Appendix  B). 

Table  10  shows  the  relative  distance  between  the  injection  location  and 
selected  gauge  locations  in  layer  6  of  the  simulation  grid. 


Table  10 

Relative  Distance  of  Gauge  Points  from  Injection  Point  for  Layer  6 

Gauge  Number 

Distance,  m 

19 

16 

32 

8 

45 

13 

100 

29 

The  observed  and  simulated  time-series  of  tritium  concentration  at 
gauge  (well)  points  19  and  32  are  shown  in  Figures  68  and  69,  respectively. 
Figures  70  and  71  compare  observed  and  simulated  concentrations  of 
naphthalene  at  gauge  points  45  and  100,  respectively. 
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Figure  68.  Comparison  of  observed  and  simulated  concentration  of  tritium  at  gauge  point  19  and  layer  6 
of  the  grid 
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Figure  69.  Comparison  of  observed  and  simulated  concentration  of  tritium  at  gauge  point  32  and  layer  6 
of  the  grid 
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Figure  70.  Comparison  of  observed  and  simulated  concentration  of  naphthalene  at  gauge  point  45  and 
layer  6  of  the  grid 
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Figure  71 .  Comparison  of  observed  and  simulated  concentration  of  naphthalene  at  gauge  point  100  and 
layer  6  of  the  grid 
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Sensitivity  Analysis 


Hydraulic  conductivity 

The  distribution  of  hydraulic  conductivity  values  has  major  effects  on 
the  flow  and  transport.  The  effect  of  decreasing  the  hydraulic  conductivity 
by  a  factor  of  10  in  the  immediate  cells  surrounding  the  injection  point  was 
investigated.  With  this  change  to  the  hydraulic  conductivities  and  the 
corresponding  leakance,  the  contaminant  plumes  demonstrated  a  tendency 
for  the  maximum  concentration  to  remain  within  a  few  meters  of  the 
injection  site.  Likewise,  the  fringes  of  the  plumes  failed  to  move  as  far 
downstream  as  before  due  to  extended  retention  time  in  the 
low-conductivity  zone. 


Porosity 

The  effect  of  changing  the  porosity  in  several  sections  of  the  aquifer  was 
evaluated.  The  MT3D  simulation  grid  was  divided  into  four  sections  for 
layers  1  through  6,  as  shown  in  Figure  72.  A  tritium  simulation  was 
conducted  with  a  three-tiered  porosity  change  to  see  the  effect  of  changing 
the  porosity.  The  porosity  in  the  area  immediately  above  the  center  line 
was  set  to  0.2  and  that  immediately  below  the  middle  line  to  0.10 
(remaining  areas  in  the  near-field  and  far-field  were  set  to  0.31).  The  effect 
of  this  was  to  move  the  front  of  the  plume  further  downgradient  than  in  the 
simulations,  with  only  a  two-tiered  or  single  value  for  porosity.  The 
simulation  conducted  with  this  three-tiered  approach  moved  the  tritium 
plume  approximately  30  m  further  downgradient  than  in  the  two-tiered 
approach.  Figure  73  shows  the  extent  of  the  tritium  plume  with  this 
three-tiered  approach  to  porosity  at  day  328. 


Recharge 

The  effect  of  recharge  rates  was  examined  in  the  MODFLOW  discussion 
(Chapter  5).  The  recharge  rate  was  increased  from  0.07  m/year  to  0.26 
m/year.  The  effect  of  the  recharge  rate  on  the  transport  results  was  minor. 
Figure  74  shows  that  the  plume’s  shape  changed  slightly  with  little  effect  on 
downgradient  extent. 


First-order  decay  rate 

The  first-order  decay  rate  (biodegradation)  was  increased  from  0.007/day 
to  0.010/day  and  MT3D/GMS  was  used  to  compute  and  illustrate 
differences  between  the  results.  The  43-percent  increase  in  the 
biodegradation  rate  caused  the  maximum  concentration  during  the  328-day 
simulation  to  decrease  by  2.5  percent.  Figure  75  shows  the  difference 
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Figure  72.  Sensitivity  analysis  using  distributed  porosities  of  0.31 ,  0.2,  0.1 ,  and 
0.31  from  top  to  bottom  in  plan  view 
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Figure  73.  Simulated  tritium  plume  distribution  at  day  328  using  porosity 
distribution  of  Figure  72 
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Figure  74.  Simulated  tritium  plume  distribution  at  day  328  with  increased 
infiltration  rate 
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Figure  75.  Naphthalene  concentration  change  due  to  an  increase  in  the  biodegradation  rate  from  0.007 
to  0.010/day  for  a)  layer  3,  b)  layer  4,  c)  layer  5,  and  d)  layer  6 
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between  a  biodegradation  rate  of  0.007/day  and  0.010/day  for  naphthalene 
for  layers  3,  4,  5,  and  6. 


Summary  and  Discussion  of  MT3D  Results 

Results  of  this  study  tend  to  indicate  that  one  of  the  most  important 
controlling  factors  in  the  simulation  of  a  heterogeneous  system  is  the 
hydrogeologic  conceptualization.  Thus,  a  realization  conditioned  on  both 
flow  and  transport  would  be  a  highly  desirable  feature  in  future  versions  of 
the  GMS.  A  better  simulation  of  the  observed  contaminant  plumes  could 
probably  be  achieved  with  a  grid  resolution  finer  than  that  used  here  (5x5 
x  1  m). 

The  simulated  transport  of  contaminants  through  the  MADE  site  was 
reasonably  close  to  the  observed  plumes,  considering  the  simplifying 
assumptions  adopted  in  the  flow  simulations.  The  flow  model  used  a 
quasi-steady-state  condition  in  which  the  piezometric  head  for  the  February 
1991  monitoring  data  set  was  used  for  the  entire  length  of  the  simulation. 
Most  of  the  simulated  conservative  tracer  (tritium)  plumes  were  close  to  the 
observed  plumes  until  approximately  234  days  into  the  simulations. 

Around  day  203,  the  average  water  table  elevation  is  known  to  increase  by 
approximately  2  m.  This  poses  significant  problems  when  attempting  to 
model  a  contaminant  plume  for  a  lengthy  period  of  time  under  these 
conditions.  The  simulated  plume  for  conservative  tracers  does  not  have  the 
required  increase  in  velocity  to  move  the  plume  nearly  the  extent  of  the 
observed  plume.  In  general,  the  tracer  monitoring  layout  for  the  site 
contained  only  a  few  measurements  above  62.5  m  which  is  the  upper  limit 
of  the  simulation  grid.  The  monitoring  layout  was  also  noted  to  increase 
the  horizontal  distance  between  measurement  locations  for  the  middle-  and 
far-field  regions.  It  is  also  noted  that  approximately  over  90  percent  of  the 
injected  contaminants  remained  above  the  injection  elevation  of  57.5  m. 
Another  observation  related  to  the  monitoring  layout  was  that  the 
monitoring  that  was  conducted  at  27  days  after  injection  used  only  a  limited 
number  of  sampling  points  near  the  injection  well.  The  observed  data 
clearly  indicate  that  in  most  of  the  observation  plots,  the  observed  plumes 
have  reached  the  furthest  extent  of  the  near-field  monitoring  extent  and 
quite  possibly  exceeded  it.  Questions  as  to  how  the  contaminants  are 
transported  in  the  region  above  an  elevation  of  62.5  m  during 
high-water-table  events  could  not  be  answered  by  MODFLOW/MT3D.  The 
limitation  of  being  able  to  model  only  the  saturated  zone  (below  the  62.5-m 
elevation)  with  MODFLOW  for  flow  and  MT3D  for  transport  is  recognized 
in  this  application  as  being  a  drawback.  In  addition,  the  lack  of  more 
detailed  contaminant  monitoring  data  in  the  upper  regions  of  the  aquifer 
above  62.5  m  may  inhibit  the  compilation  of  an  exhaustive  data  set  for  the 
MADE  site  that  can  be  used  as  a  test  data  set  for  the  complete  analysis  of 
present  simulation  models.  The  MADE  data  set  is  an  impressive  one  and  is 
the  most  complete  for  a  heterogeneous  site  to  date. 
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Previous  reviews  of  the  MADE  experiments  (Boggs  et  al.  1993) 
indicated  that  there  is  a  strong  possibility  that  organic  tracer  concentrations 
could  have  been  significantly  reduced  by  the  sampling  mechanics  due  to 
organic  tracer  volatility.  Most  of  the  organic  tracer  simulations  indicated  a 
plume  that  traveled  further  downgradient  than  the  observed  plume  toward 
the  latter  half  of  the  simulations.  The  exception  to  this  was  the 
o-dichlorobenzene  simulation,  in  which  the  observed  plume  traveled  further 
downgradient  than  the  simulated  plume  and  showed  extensive  spatial 
discontinuity.  The  p-xylene  simulation  demonstrated  the  transformation  of 
the  chemical  as  compared  to  the  14C  plume  which  traveled  further 
downgradient  than  the  organic  chemical.  This  is  noted  because  the  p-xylene 
was  labeled  with  the  14C  that  was  tracked  in  the  monitoring. 

The  simulations  were  very  sensitive  to  hydraulic  conductivity  values.  A 
change  in  the  hydraulic  conductivity  at  the  injection  well  caused  major 
effects  for  the  travel  of  the  contaminants,  as  expected. 

The  simulation  did  not  appear  to  be  very  sensitive  to  the  recharge  rates 
evaluated  in  the  sensitivity  analysis.  The  increase  in  the  recharge  rate  in  the 
flow  model  from  0.07  m/year  to  0.26  m/year  induced  only  a  small  increase 
in  the  tritium  plume  and  the  shape  of  that  plume. 

The  model  was  sensitive  to  changes  in  porosity.  A  three-tiered  porosity 
profile  that  followed  the  meander  channel  (with  a  porosity  value  in  the 
far-field  region  of  0.31,  the  middle-field  region  of  0.2,  in  the  meander 
channel  of  0.1,  and  in  the  near-field  region  of  0.31)  caused  the  conservative 
tracer  tritium  plume  to  move  approximately  30  m  further  than  the  tritium 
plume  simulated  with  a  porosity  profile  of  0.10  in  the  meander  channel  and 
0.31  elsewhere. 

Recommendations  for  future  simulations  include  better  representations 
of  the  hydraulic  conductivity  for  the  site  and  a  model  grid  with  finer 
resolution. 


116 


Chapter  6  MT3D  Modeling  Approach 


References 


Boggs,  J.  M.,  and  Young,  S.C.  (1988).  “Evaluation  of  tracer  sampling 
devices  for  macrodispersion  experiment,”  EPRI  Report  EA-5816, 
Electric  Power  Research  Institute,  Palo  Alto,  CA. 

Boggs,  J.  M.,  Beard,  L.  M.,  and  Waldrop,  W.  R.  (1993).  “Transport  of 
tritium  and  four  organic  compounds  during  a  natural-gradient 
experiment  (MADE-2),”  EPRI  Report  TR-101998,  Electric  Power 
Research  Institute,  Palo  Alto,  CA. 

Boggs,  J.  M.,  Beard,  L.  M.,  Waldrop,  W.  R.,  Stauffer,  T.  B.,  MacIntyre, 
W.  G.,  and  Antworth,  C.  P.  (1993).  “Transport  of  tritium  and  four 
organic  compounds  during  a  natural-gradient  experiment  (MADE-2),” 
Technical  Report,  EPRI  TR-101998,  Project  2485-05,  Palo  Alto,  CA. 

Boggs,  J.  M.,  Young,  S.  C.,  Benton,  D.  J.,  and  Chung,  Y.  C.  (1990). 
“Hydrogeologic  characterization  of  the  MADE  site,”  EPRI  Report 
EN-6915,  Electric  Power  Research  Institute,  Palo  Alto,  CA. 

Boggs,  J.  M.,  Young,  S.  C.,  Hopkins,  R.  A.,  and  Chung,  Y.  C.  (1988). 
“Interim  hydrogeologic  characterization  of  the  Macrodispersion 
Experiment  Site,”  Draft  Technical  Report  for  EPRI,  Research  Project 
2485-5,  Palo  Alto,  CA. 

Brigham  Young  University.  (1995).  “The  Department  of  Defense 
Groundwater  Modeling  System,”  Engineering  Computer  Graphics 
Laboratory,  Brigham  Young  University,  Provo,  UT. 

Electric  Power  Research  Institute.  (1989).  “Estimation  of 

macrodispersivity  based  on  analysis  of  hydraulic  conductivity 
variability  at  the  MADE  site,”  EPRI  EN-6405,  Project  2485-5,  Palo 
Alto,  CA. 


_ .  (1990).  “Hydrogeologic  characterization  of  the  MADE  site,” 

EPRI  EN-6915,  Project  2485-5,  Palo  Alto,  CA. 

_ .  (1991).  “Database  for  the  first  macrodispersion  experiment 

(MADE-1),”  EPRI  EN-7363,  Project  2485-5,  Palo  Alto,  CA. 


References 


117 


Electric  Power  Research  Institute.  (1993a).  “Database  for  the  second 
macrodispersion  experiment  (MADE-2),”  EPRI  TR-102072,  Project 
2485-05,  Palo  Alto,  CA. 

_ .  (1993b).  “Transport  of  tritium  and  four  organic  compounds 

during  a  natural-gradient  experiment  (MADE-2),”  EPRI  TR-101998, 
Project  2485-05,  Palo  Alto,  CA. 

Gray,  D.  D.,  and  Dale,  F.  R.  (1995).  “Transport  modeling  of  tritium  at  the 
MADE-2  site.”  Proceedings  of  Groundwater  Management.  R.  J. 
Charbeneau,  ed.,  American  Society  of  Civil  Engineers,  New  York, 
164-168. 

Hill,  M.  C.  (1990).  “MODFLOWP:  A  computer  program  for  estimating 
parameters  of  a  transient,  three-dimensional  groundwater  flow  model 
using  nonlinear  regression,”  Open-File  Report,  U.S.  Geological  Survey. 

Leake,  S.  A.,  and  Prudic,  D.  E.  (1988).  “Documentation  of  a  computer 
program  to  simulate  aquifer-system  compaction  using  the  modular 
finite-difference  ground-water  flow  model,”  Open-File  Report  88-482, 
Appendix  C,  69-80,  U.S.  Geological  Survey,  AZ. 


Lin,  H-C.  J.,  Richards,  D.  R.,  Yeh,  G-T.,  Cheng,  J-R.,  Cheng,  H-P.,  and 
Jones,  N.  L.  (1997).  “FEMWATER:  A  three-dimensional  finite  element 
computer  model  for  simulating  density-dependent  flow  and  transport  in 
variably  saturated  media,”  Technical  Report  CHL-97-12,  U.S.  Army 
Engineer  Waterways  Experiment  Station,  Vicksburg,  MS. 

McDonald,  M.  G.,  and  Harbaugh,  A.  W.  (1988).  “A  modular 

three-dimensional  finite-difference  ground-water  flow  model,”  U.S. 
Geological  Survey  Techniques  of  Water-Resources  Investigations, 

Book  6,  Chapter  Al. 

Mullins,  J.  A.,  Carsel,  R.  F.,  Scarbrough,  J.  E.,  and  Ivery,  A.  M.  (1993). 
“PRZM-2,  a  model  for  predicting  pesticide  fate  in  the  crop  root  and 
unsaturated  soil  zones:  Users  manual  for  release  2.0,”  Environmental 
Research  Laboratory,  U.  S.  Environmental  Protection  Agency,  Athens, 
GA. 

Rehfeldt,  K.  R.,  Boggs,  J.  M.,  and  Gelhar,  L.  W.  (1992).  “Field  study  of 
dispersion  in  a  heterogeneous  aquifer;  3,  Geostatistical  analysis  of 
hydraulic  conductivity,”  Water  Resources  Research. 

Rehfeldt,  K.  R.,  Gelhar,  L.  W.,  Southard,  J.  B.,  and  Dasinger,  A.  M. 
(1989b).  “Estimates  of  macrodispersivity  based  on  analysis  of 
hydraulic  conductivity  variability  at  the  MADE  site,”  EPRI  Report 
EN-6405,  Electric  Power  Research  Institute,  Palo  Alto,  CA. 


118 


References 


Rehfeldt,  K.  R.,  Hufschmied,  P.,  Gelhar,  L.  W.,  and  Schaefer,  M.  E. 
(1989a).  “The  borehole  flowmeter  technique  for  measuring  the  vertical 
distribution  of  hydraulic  conductivity,”  Electric  Power  Research 
Institute,  Technical  Report  EN-6511,  Palo  Alto,  CA. 

Zheng,  C.  (1990).  “MT3D:  A  modular  three-dimensional  transport  model 
for  simulation  of  advection,  dispersion  and  chemical  reactions  of 
contaminants  in  groundwater  systems,”  S.S.  Papadopulos  and 
Associates,  Rockville,  MD. 


References 


119 


Appendix  A 

Interpolated  Hydraulic  Conductivity 


The  interpolated  hydraulic  conductivity  values  (units  of  m/day)  for  each  layer  of  the  model 
grid  are  shown  in  the  following  Figures  Al-All. 
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Figure  A1 .  Interpolated  Hydraulic  Conductivity  for  Model  Layer  1 
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Figure  A2.  Interpolated  Hydraulic  Conductivity  for  Model  Layer  2 
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Figure  A3.  Interpolated  Hydraulic  Conductivity  for  Model  Layer  3 
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Figure  A4.  Interpolated  Hydraulic  Conductivity  for  Model  Layer  4 
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Figure  A5.  Interpolated  Hydraulic  Conductivity  for  Model  Layer  5 
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Figure  A6.  Interpolated  Hydraulic  Conductivity  for  Model  Layer  6 
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Figure  A7.  Interpolated  Hydraulic  Conductivity  for  Model  Layer  7 
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Figure  A8.  Interpolated  Hydraulic  Conductivity  for  Model  Layer  8 
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Figure  A9.  Interpolated  Hydraulic  Conductivity  for  Model  Layer  9 
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Figure  A10.  Interpolated  Hydraulic  Conductivity  for  Model  Layer  10 
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Figure  A1 1 .  Interpolated  Hydraulic  Conductivity  for  Model  Layer  1 1 
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Appendix  B 

Piezometric  Head  Contours 


This  appendix  shows  the  average  head  contours  for  both  the  shallow  and  deep  parts  of  the 
aquifer  for  each  day  of  the  available  data.  The  plots  follow  and  are  not  labeled  with  individual 
figures  since  there  are  several  and  the  headings  are  believed  to  be  self-explanatory. 
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Appendix  C 

Gauge  Plots  of  Transient  Flow  Simulations 


This  appendix  shows  some  of  the  Groundwater  Modeling  system  (GMS)  gauge  plot 
comparisons  for  an  attempt  at  a  transient  flow  simulation  with  a  grid  consisting  of  over  175,000 
cells.  The  plots  show  the  correlation  between  the  simulated  head  and  the  observed  head  for  the 
simulation.  Transient  simulations  using  this  enhanced  grid  required  over  1  gigabyte  of  storage 
for  a  single  output  file.  The  attempt  was  aborted  due  to  a  combination  of  problems  including 
memory,  space,  unsaturated  zones  in  the  model,  too  many  boundary  conditions  driving  the 
simulation,  etc. 
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improved  overall  by  simulating  a  low-  and  high-water  table  for  the  simulation  period  of  328  days  as  described  in  the 
FEMWATER  application  chapters.  The  simulation  without  the  high-water-table  event  showed  that  the  travel  extent 
of  simulated  contaminant  plumes  is  limited  compared  to  observed  plumes.  Most  of  the  MT3D  simulated 
conservative  tracer  (tritium)  plumes  were  close  to  the  observed  plumes  until  approximately  234  days  into  the 
simulations.  Previous  reviews  of  the  MADE  experiments  indicated  that  there  is  a  strong  possibility  that  the  organic 
tracer  concentrations  could  have  been  significantly  reduced  by  volatility  losses.  Most  of  the  organic  tracer 
simulations  indicated  a  plume  that  traveled  further  downgradient  than  the  observed  plume  toward  the  latter  half  of 
the  simulations.  The  exception  to  this  was  the  o-dichlorobenzene  simulation  in  which  the  observed  plume  traveled 
further  downgradient  than  the  simulated  plume  and  showed  extensive  spatial  discontinuity. 


